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PREFACE 


This  report  includes  the  results  of  two  studies  with  implications  for 
yield  estimation  of  underground  nuclear  explosions.  In  the  first  part,  we 
develop  a  method  for  the  simultaneous  inversion  of  near-field  waveforms  for 
source  and  structure  parameters.  This  is  desirable  in  order  to  determine  the 
extent  to  which  assumptions  of  the  velocity  structure  affect  estimates  of 
source  parameters.  The  simultaneous  inversion  automatically  provides 
quantitative  information  on  parameter  trade-offs,  as  well  as  error  analyses  on 
the  resultant  ncdels.  The  inversion  also  provides  a  means  of  determining 
source  parameters  at  sites  where  no  structure  model  previously  exists,  as 
would  be  the  case  in  an  on-site  verification  experiment.  Test  inversions  are 
presented,  as  well  as  preliminary  results  of  inversions  of  vertical  velocity 
waveforms  from  Pahute  Mesa  events.  The  second  part  of  this  report  addresses 
the  implications  of  choosing  a  particular  representation  of  the  RDP  on  the 
determination  of  the  teleseisraic  attenuation  operator.  Using  source 
parameters  determined  from  near-field  modeling  assuming  either  the 
Mueller-Murphy  or  Helmberger-Hadley  effective  source  functions,  time-domain 
amplitude  and  period  information  from  a  large  set  of  short-period,  WWSSN 
teleseismic  recordings  were  modeled  to  determine  both  frequency-independent 
and  frequency-dependent  Q  models.  Including  spectral  shape  data  in  the  0.5  to 
4  Hz  band  as  constraints,  either  source  model,  when  convolved  with  a 
frequency-dependent  Q  operator,  is  consistent  with  both  near-field  and 
teleseismic  observations. 


ABSTRACT 


The  determination  of  source  parameters  for  underground  nuclear  explosions 
through  waveform  modeling  of  near-field  recordings  is  dependent  upon  the 
velocity  structure  model  assumed  between  the  source  and  receivers.  In  order 
to  minimize  errors  due  to  inadequacies  in  the  assumed  structure  model,  we  have 
developed  a  simultaneous  inversion  of  near-field  waveforms  for  source  and 
structure  parameters.  The  inversion  automatically  provides  information  on  the 
trade-offs  between  parameters  and  error  analyses  of  the  resultant  models.  The 
source  may  be  parameterized  by  any  of  several  effective  source  functions, 
although  for  this  study  we  use  the  representation  of  Helmberger  and  Hadley 
(I98l).  To  minimize  free  parameters,  the  structure  is  parameterized  as  a 
layered  stack  of  gradients  defined  by  the  velocity  at  the  top  of  the  layer, 
the  velocity  gradient  within  the  layer  and  the  depth  to  the  top  of  the  layer. 
The  gradient  model  is  then  discretized  into  homogeneous  plane  layers,  with  the 
layer  thickness  defined  as  the  wavelength  of  the  predominant  period  in  the 
data,  and  synthetics  are  computed  using  asymptotic  generalized  ray  theory. 

For  each  observed  waveform,  the  residuals  consist  of  the  normalized 
cross-correlation  coefficient  between  the  data  and  synthetics  (a  measure  of 
waveform  fit),  the  difference  in  the  normalization  factors  (a  measure  of 
absolute  amplitude),  and  the  time  lag  to  the  peak  of  cross-correlation  (a 
measure  of  absolute  travel  time).  Numerical  partial  derivatives  are  computed 
by  perturbing  the  starting  model,  and  the  problem  is  solved  using  an  iterative 
generalized  inverse. 

Four  test  inversions  of  synthetic  data  sets  are  presented  which 
illustrate  the  need  to  truncate  eigenvalues  and  preferentially  weight  the 


parameters  and  residuals  in  order  to  stabilize  the  inversion.  Generally  in 
these  tests,  when  the  parameters  are  well  resolved,  the  inversion  is  quite 
robust  and  convergence  is  rapid.  Preliminary  results  are  presented  for 
inversions  of  the  vertical  velocity  waveforms  from  the  Pahute  Mesa  events, 
BOXCAR,  INLET,  MAST  and  SCOTCH.  For  this  report,  the  structure  models  are 
constrained  to  two  gradients  separated  by  a  second-order  discontinuity.  The 
structure  model  obtained  for  BOXCAR  is  not  significantly  different  than  the 
Pahute  Mesa  model  of  Hartzell,  et  al.  (1983),  but  the  models  obtained  for 
INLET,  MAST,  and  SCOTCH  suggest  a  higher  gradient  in  the  upper  layer,  and  poor 
resolution  of  the  lower  layer.  The  source  models  obtained  in  these 
preliminary  inversions  give  higher  estimates  of  V  than  those  of  Barker,  et 
al.  (1985). 


INTRODUCTION 


The  determination  of  source  parameters  for  underground  nuclear  explosions 
through  waveform  modeling  of  near-field  recordings  is  dependent  upon  the 
velocity  structure  model  assumed  between  the  source  and  receivers.  While 
independent  information  may  exist  for  some  aspects  of  the  structure,  such  as 
the  velocity  at  the  surface  and  the  depths  to  the  water  table  and  basement, 
the  detailed  models  are  generally  the  result  of  laborious  tr ial-and-error 
modeling  of  observed  waveforms  (e.g.,  Helmberger  and  Hadley,  1981;  Hartzell, 
et  al.  ,  1983;  Burdick,  et  al. ,  1984).  Discretizations  of  velocity  gradients 
into  plane-layered  structure  models  must  be  fine  enough  to  avoid  introducing 
high  frequency  errors  into  the  source  function  determinations.  For  signals 
with  frequencies  up  to  5  Hz  and  compressional  velocities  near  3  km/s,  layer 
thickness  should  be  less  than  about  0.6  km.  This  is  particularly  important 
near  the  "'■'urce  depth  due  to  the  interference  of  upgoing  and  diving  energy 
over  near-field  ranges.  Thus  an  adequate  structure  model  for  near-field 
waveform  modeling  might  consist  of  more  than  ten  layers  in  the  top  10  km  of 
the  crust.  There  is  obviously  a  premium  to  be  paid  before  undertaking  the 
modeling  of  sources  at  different  test  sites.  Further,  the  structure  model 
obtained  is,  to  a  certain  extent,  dependent  upon  the  source  function  used 
during  the  modeling  procedure.  The  result  is  clearly  non-unique  and  the 
extent  of  the  trade-offs  cannot  easily  be  quantified. 

One  solution  to  this  problem  is  a  simultaneous  inversion  for  source  and 
structure  parameters.  The  complexity  of  wave  interaction  in  the  near-field 
makes  such  an  inversion  an  approximate  and  somewhat  costly  endeavor. 
Nevertheless,  while  the  models  obtained  by  such  an  inversion  may  not  be 
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intrinsically  "better"  than  their  tr ial-and-error  counterparts,  a  suitable 
choice  of  inversion  method  will  automatically  proviie  invaluable  information 
concerning  parameter  and  data  resolution.  Not  only  can  the  source  vs. 
structure  trade-offs  be  quantified  in  this  way,  but  also  the  trade-offs  within 
the  source  and  structure  parameterizations  (such  as  the  B-  trade-off  in  the 
RDP)  may  be  determined.  We  also  obtain  an  indication  of  which  aspects  of  the 
observed  data  are  most  important  in  constraining  the  solution. 

This  section  of  this  report  includes  discussions  of  the  development  of 
such  an  inversion,  test  inversions  using  synthetic  data  sets,  and  applications 
to  near-field  data  from  the  Pahute  Mesa  test  site.  Following  reports  will 
include  inversions  of  data  from  other  test  sites,  as  well  as  modifications  of 
the  assumed  parameterizations  for  the  Pahute  Mesa  events. 


BACKGROUND 


The  vertical  displacement  due  to  an  explosion  source  may  be  written  using 
first  order  asymptotic  generalized  ray  theory  as  (see  for  example,  Helmberger 
and  Harkrider,  19TB) 


w(r,t)  =  4*  (t) 
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where  *  denotes  convolution,  4,,(t)  is  the  time  derivative  of  the  reduced 
displacement  potential  (RDP),  and 


R.(p)  =  ^  ni(P>  dt  ^  c 

which  is  evaluated  along  the  Cagniard  contour,  p  is  the  ray  parameter,  na  is 
the  vertical  wave  slowness  ( n  =  /a  -  p^  ),  R^z(p)  is  the  receiver  function 
defined  by  Helmberger  (1971*),  and  TI^p)  is  the  product  of  reflection  and 
transmission  coefficients  along  the  path  for  each  ray. 

For  the  inversion,  we  need  to  compute  the  partial  derivatives  of  (l)  with 
respect  to  the  individual  parameters  of  the  source  and  velocity  structure. 

The  partial  with  respect  to  structure  parameters,  say  the  velocity  of  the 
layer,  Vj ,  reduces  to  the  evaluation  of 

9R.  (p)  ^(P)  9P  (3) 

~3^T  “  3p  3v 

3  3 

which  consists  of  evaluating  the  partials  of  the  various  terms  in  (2)  as  a 
function  of  the  complex  ray  parameter  along  the  Cagniard  contour,  while  the 
contour  itself  varies  with  changes  in  the  velocity  parameter.  This  evaluation 
must  be  made  for  each  ray,  then  summed  to  obtain  the  partial  derivative. 

Previous  structure  inversions  based  on  waveform  modeling  have  invoked  a 
variety  of  approximations  in  order  to  simplify  (3).  Mellman  (19B0)  utilized  a 
modified  first  motion  approximation,  such  that  rapidly  varying  terms  in  (2) 
are  approximated  by  evaluating  them  only  at  the  geometric  ray  parameter.  Shaw 
(19^3)  and  Chapman  and  Orcutt  (1985)  used  the  WKBJ  approximation  to  invert 
marine  reflection/refraction  data,  and  Brown  (19S2)  employed  a  similar 
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approach  to  obtain  the  ocean  sound  speed  profile.  Since  the  WKBJ 
approximation  becomes  invalid  where  large  velocity  gradients  exist,  or  where 
discontinuities  are  encountered  near  the  turning  depth  of  a  ray.  Given  ( I98U ) 
developed  a  hybrid  WKBJ-generalized  ray  approach  to  obtain  upper  mantle 
structure.  In  inverting  long-period,  regional  Pn^  waveforms  for  average 
crustal  structure,  Wallace  (1983)  generated  average  generalized  ray  Green's 
functions,  which  were  "stretched"  or  "squeezed"  to  simulate  perturbations  in 
the  structure  model. 

In  forward  modeling  of  near-field  waveforms,  Burdick,  et  al.  (1982)  and 
Barker,  et  al.  (1985)  have  shown  that  the  observed  waveforms  result  from  the 
interference  of  upgoing  and  diving  rays,  as  well  as  surface  reflections.  This 
interference  is  particularly  sensitive  to  the  velocity  gradient  near  the 
source  depth.  Each  synthetic  seismogram  consists  of  a  superposition  of  rays 
from  a  variety  of  take-off  angles  that  sample  different  portions  of  the 
structure.  An  algorithm  designed  to  determine  if  and  when  the  modified  first 
motion,  WKB.T  and  other  approximations  are  appropriate  for  a  particular  ray,  to 
compute  the  analytic  partial  derivatives,  and  sum  over  rays  would  be  of 
questionable  efficiency  and  accuracy.  Instead,  numerical  partials  may  be 
computed  by  evaluating  the  terms  in  (3)  along  the  complex  ray  parameter  for 
each  ray.  This  approach  is  obviously  nonlinear,  and  care  must  be  taken  to 
stabilize  the  resulting  parameter  changes.  For  simpler  problems,  the 
computation  of  numerical  partial  derivatives  would  be  quite  inefficient,  but 
given  the  complexity  of  the  near-field  waveform  modeling  problem,  the  accuracy 
of  the  method  will  more  than  counterbalance  the  computational  cost. 

A  number  of  techniques  have  been  used  to  obtain  source  parameters  by 
inversion  of  body  waveforms.  Burdick  and  Mellman  (1976)  inverted  teleseismic 
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body  waves  to  obtain  the  orientation  and  time  history  of  the  Borrego  Mt. 
earthquake.  Several  investigators  have  inverted  teleseismic  body  waves  to 
obtain  the  seismic  moment  tensor.  A  few  examples  include  Stump  and  Johnson 
(1977),  Langston  (1981),  Dziewonski,  et  al.  (1981)  and  Nabelek  (1984).  To  our 
knowledge,  the  first  inversion  of  near-field  data  for  explosion  source 
parameters  was  performed  by  Stump  and  Johnson  (1984).  McEvilly  and  Johnson 
(1986)  and  Wallace  ( 1986 )  have  recently  reported  on  similar  applications.  In 
these  studies,  the  six  independent  elements  of  the  moment  tensor  are  obtained 
and  the  relative  amplitudes  and  time  histories  of  the  isotropic  and  deviatoric 
parts  are  interpreted  in  terms  of  explosion  and  tectonic  release  sources, 
respectively.  Because  the  moment  tensor  is  free  to  take  on  any  value,  and 
because  the  velocity  structure  model  is  fixed,  it  is  entirely  possible  that 
errors  in  the  structure  model  and  noise  in  the  data  could  map  into  the 
non-isotropic  moment  tensor.  In  addition,  since  no  constraints  are  placed  on 
the  form  of  the  isotropic  time  history,  it  is  not  easily  interpreted  in  terms 
of  the  parameters  of  common  representations  of  the  explosion  RDP. 

Because  of  our  assertion  that  errors  in  the  velocity  structure  adversely 
affect  explosion  source  parameter  estimates,  we  are  interested  in  approaching 
these  problems  before  tackling  source  complexity.  Therefore,  our  preliminary 
inversions  will  constrain  the  explosion  source  to  be  isotropic,  and  we  will 
constrain  the  source  time  history  to  have  the  form  of  the  time  derivative  of 
the  RDP.  Various  forms  of  the  RDP  may  be  considered,  and  the  ramifications  of 
that  choice  will  be  investigated  in  future  work.  After  evaluating  these 
results,  the  source  parameterization  could  be  generalized  to  allow  the  source 
time  history  to  take  on  an  arbitrary  shape.  At  a  later  stage,  a  deviatoric 
moment  tensor  source  representing  tectonic  release  could  be  added  to  the 
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isotropic  source,  with  the  latter  constrained  to  the  known  depth  and  origin 
time  of  the  explosion.  Inverting  for  a  single,  six-component  moment  tensor, 
however,  constrains  the  non-isotropic  portion  of  the  source  to  occur  at  the 
same  location  as  the  explosion.  While  this  may  turn  out  to  be  the  case, 
near-field  observations  should  allow  resolution  of  spatial  separations  in  the 
two  sources,  so  such  a  constraint  may  be  unreasonable. 


INVERSION  METHOD 


The  layered  structure  models  required  for  generalized  ray  synthetics  may 
be  considered  discretizations  of  one  or  more  velocity  gradients  separated  by 
first-  or  second-order  discontinuities.  In  order  to  minimize  the  number  of 
free  parameters  in  the  Inversion,  we  will  utilize  this  fact  and  solve  not  for 
the  velocities  and  thicknesses  of  individual  layers,  but  for  the  velocity 
gradient  and  the  velocity  at  the  top  of  the  gradient.  As  shown  In  Figure  1, 
the  discretized  model  follows  this  gradient,  with  layer  thicknesses  determined 
by  the  wavelength  of  the  predominant  period  in  the  data.  Shear  velocity  and 
density  may  be  free  parameters,  or  may  be  related  to  the  corapressional 
velocity  by  predefined  constants.  In  general,  the  free  parameters  are  the 
gradient  within  each  layer,  the  velocity  at  the  top  of  the  layer,  and  the 
depth  to  the  top  of  the  layer.  Thus,  for  an  n-layered  structure  in  which 
shear  velocity  and  density  are  not  free  parameters,  the  number  of  structure 
parameters  in  the  inversion  is  3n-l.  The  number  of  gradients  to  be  considered 


will  be  fixed  for  each  Inversion,  so  the  optimal  number  will  have  to  be 

'"’mined  by  test  inversions  proceeding  from  simpler  to  more  complex  models 


DEPTH  (km) 


Structure  Parameterization 


Figure  1.  Structure  model  parameterization.  The  layered  structure  model  is 
discretized  from  one  or  more  velocity  gradients  defined  by  the  values  of 
the  gradient,  the  velocity  at  the  top  of  the  gradient  and  the  depth  to 
the  top  of  the  gradient.  The  depth  increment  is  defined  as  the  product 
of  the  P-vave  velocity  and  the  predominant  period  of  the  data  (typically 
0.25  sec).  S-wave  velocity  and  density  are  linearly  related  to  the 
P-vave  velocity.  The  case  shovn  consists  of  two  gradients  separated  by  a 
second-order  discontinuity. 
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Constraints  such  as  the  depth  to  the  water  table,  or  that  a  particular 
interface  nay  have  only  a  second-order  discontinuity,  reduce  the  number  of 
free  parameters  while  increasing  non-linearity  of  the  inversion. 

Since  the  discretized  velocity  model  will  vary  from  iteration  to 
iteration  as  well  as  with  the  model  perturbations  required  to  compute  the 
numerical  partials,  the  ray  set  from  which  the  synthetics  are  computed  will 
also  vary.  Barker,  et  al.  (1985)  found  that  the  first  one  or  two  cycles  of 
vertical  velocity  observations  at  near-field  ranges  were  well  modeled  by 
including  the  upgoing  direct  ray,  a  sum  of  diving  rays  that  constitute 
turning,  reflecting  and  critically  refracting  energy ,  and  a  sum  of  rays  that 
depart  upward,  reflect  from  the  free  surface  and  follow  diving  ray  paths.  If 
radial  observations  are  to  be  modeled,  or  if  later  portions  of  the  waveform 
are  considered,  P-to-S  conversions  at  the  free  surface  and  at  discontinuities 
should  be  included.  An  algorithm  to  generate  this  ray  set  for  each  new 
layered  structure  is  simple  to  implement,  and  computational  time  is  saved  by 
culling  rays  which  arrive  outside  the  inversion  time  window,  or  whose 
amplitude  will  be  below  some  predefined  cutoff  level. 

For  this  report,  the  source  parameterization  follows  the  form  of  the  RDP 
proposed  by  Helmberger  and  Hadley  (1981): 

^(t)  =  ^{l  -  e"Kt  [1  +  Kt  +  J( Kt)2  -  B(Kt)U  (4) 

where  ^represents  the  static  value,  K  represents  the  rise  time,  and  B 
represents  the  overshoot.  Haskell  (1967)  proposed  that  the  polynomial  in  (U) 
should  be  quartic,  while  Mueller  and  Murphy  (1971)  and  von  Seggern  and 
Blandford  (1972)  suggest  quadratic  forms.  Barker,  et  al.  (1Q86)  showed  that  a 
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cubic  RDP  convolved  with  elastic  Green's  functions  explains  near-field 
observations.  A  quadratic  RDP  can  explain  the  near-field  observations  equally 
well  as  long  as  the  Green's  functions  include  the  effects  of  depth-dependent 
crustal  attenuation.  Although  attenuation  could  be  included  in  the  inversion, 
this  would  greatly  increase  the  number  of  free  parameters,  and  attenuation 
parameters  would  likely  be  poorly  resolved.  Initially  therefore,  we  choose 
the  Helmberger-Hadley  RDP  and  elastic  wave  propagation  as  the  simpler  model. 
The  effective  source  function  may  be  obtained  by  taking  the  first 
(displacement),  second  (velocity)  or  third  (acceleration)  time  derivatives  of 
the  RDP.  As  may  be  seen  from  (U),  the  RDP  and  its  time  derivatives  are  linear 
-n  and  B,  but  not  in  K.  Burdick,  et  al.  (1982)  indicated  that  4^  and  B 
trade  off  nearly  totally  in  modeling  near-field  velocity  records.  While  this 
trade-off  may  be  quantified  by  allowing  both  parameters  to  vary  in  the 
inversion,  useful  source  estimates  may  require  that  B  be  constrained. 

The  linearized  least-squares  inverse  may  be  obtained  by  solving 

Ac'  =  A'  Ap'  (5) 

in  which  Ac'  is  the  residual  vector,  containing  the  differences  between 
observed  and  synthetic  data.  A'  is  the  matrix  of  partial  derivatives,  and  Ap' 
is  the  desired  vector  of  parameter  changes.  Perhaps  the  simplest 
determination  of  the  residual  vector  is  a  point-by-point  difference  of 
observed  and  synthetic  seismograms  within  the  time  window.  This  implies, 
however,  that  neighboring  points  in  a  seismogram  represent  independent 
observations,  which  is  generally  not  the  case.  On  the  other  hand,  since  each 
seismogram  reflects  a  lagged  sum  of  rays  which  sample  different  portions  of 
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the  structure,  the  resolution  of  parameters  varies  through  the  time  window, 
and  a  single  observation  for  each  seismogram  underestimates  the  information 
contained  in  the  waveform. 

Burdick  and  Mellraan  (1976)  utilized  an  error,  or  residual,  function  based 
on  the  normalized  cross-correlation  coefficient  between  observed  and  synthetic 
waveforms.  That  is, 

ej  =  1  -  max(o^  cc  s^)  (6) 


where 
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o .  (t) 
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o.  (t) 
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(Jo*<t>  dt)’-/2 


(7) 


and 


s.  (t) 
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s.  (t) 


(Js^(t)  dt)l/z 


(8) 


o^(t)  is  the  ith  observed  seismogram,  s^(t)  is  the  ith  synthetic  seismogram 
and  cc  denotes  cross-correlation.  This  error  function  has  become  popular 
because  it  is  sensitive  to  the  shapes  of  the  waveforms,  while  the 
normalization  makes  it  insensitive  to  absolute  amplitude,  and  measuring  the 
maximum  of  the  cross-correlation  makes  it  insensitive  to  absolute  time.  In 
the  near-field  problem,  however,  absolute  amplitude  and  time  are  well 


calibrated,  and  should  be  fit  by  the  resultant  model.  Therefore,  in  addition 


> 
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to  (6),  our  residual  will  include  for  each  seismogram  the  relative  difference 
in  the  normalization  factors. 


(Jju)  dt)Vi  -  (/,*(t>  dt),/! 

dt)  /! 


which  is  a  measure  of  absolute  amplitude  residual,  and  the  time  lag  to  the 
maximum  of  the  cross-correlation,  tj,  which  is  a  measure  of  absolute  time 
residual.  The  residual  vector  is  now  defined  as 


Ac'  = 


Sl*  Rl*  tl*  a2’  t2 . em*  am»  tj' 


for  m  seismograms  and  the  objective  function  to  be  minimized  by  the  inverse  is 

r  -  I  <e;2  .  .  tf)  .  (li) 

m 

With  this  choice  for  the  residual  function,  the  numerical  partial  derivatives 
become,  rather  than  (3), 

8ei  (ei<Pi+Api)  -  e.(p.-Ap.)) 

3Pd  ~  2  Ap.  (12) 

1  J 

where  APj  is  the  perturbation  to  the  starting  model  parameter,  Pj.  The 
partials  ;*ai/JPj  and  ■»ti/^pJ  have  similar  form. 

Now  that  we  have  defined  the  form  of  the  residual  and  parameter  change 
vectors,  we  may  redefine  the  mtrices  in  (5)  as 
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Ac  =  s-1^2  Ac' 


Ap  =  v-^/2  Ap* 


A  =  s -1'2  A  w1/2 


in  which  s  is  the  covariance  matrix  of  the  residuals  and  w  is  a  weighting 
matrix  of  the  parameter  changes.  For  simplicity,  we  assume  that  errors  are 
uncorrelated,  so  s  is  diagonal.  If  an  estimate  of  the  residual  variance  is 
known  (e.g.  ,  timing  errors),  that  value  is  used  in  s.  Otherwise,  the  residual 
variances  are  inversely  related  to  the  relative  confidence  in  the  residuals. 
The  basic  purpose  of  the  diagonal  matrix  w  is  to  counter  the  dimensionality  of 
the  parameterization  of  source  and  structure.  Obviously  velocity  gradient 
(with  values  on  the  order  of  1  km/s/km)  and  ^  (values  of  1010  cm^)  should  not 
be  weighted  equally.  Due  to  nonlinearity,  the  weighting  factors  are  also  used 
to  stabilize  the  parameter  changes.  For  example,  if  a  parameter  changes  to  an 
unreasonable  value  or  oscillates  about  some  value,  the  weight  of  that 
parameter  is  decreased  and  the  inversion  is  repeated.  This  introduces  a 
certain  degree  of  uncertainty  in  interpreting  the  variances  of  the  parameter 
changes,  since  these  now  include  the  effects  of  the  somewhat  arbitrary  weights 
necessary  due  to  nonlinearity.  Although  convergence  criteria  may  be  utilized 
to  determine  when  to  stop  an  inversion,  we  have  not  implemented  such  an 
approach,  but  simply  allow  additional  sets  of  five  iterations  until 
convergence  is  deemed  complete  or  unattainable.  Finally,  following  Jackson 
(1979) »  a  priori  constraints  may  be  placed  on  the  parameter  changes  so  that, 
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for  example,  the  velocity  cannot  be  negative,  or  90  that  source  parameters 
must  fall  within  a  range  of  reasonable  values. 

The  weighted  version  of  (5)  may  be  solved  by  a  linearized  generalized 
inverse  (Wiggins,  1972)  such  that 

Ap  =  A+  Ac  =  V  A-1  UT  Ac  (14) 

in  which  A  is  the  diagonal  matrix  of  non-zero  singular  values  of  A,  V  is  the 
corresponding  matrix  of  eigenvectors  which  spans  the  parameter  space,  and  U  is 
the  matrix  of  eigenvectors  that  spans  the  residual  space.  Since  the  problem 
is  non-linear,  the  inversion  is  Iterative  and  minimizes  both  the  least-squares 
error,  |a  Ap  -  Ac  p,  and  the  length  of  the  parameter  changes,  [ Ap |  .  Since 
small  singular  values  result  in  large  parameter  changes,  an  a  priori  cut-off 
is  specified,  below  which  the  singular  value  is  assumed  to  be  zero,  and  the 
dimension  of  the  parameter  space  is  truncated.  The  variances  of  the  parameter 
changes  are  given  by 

O?  -  I  vfj  /  ^ 

j 

and  include  the  effects  of  the  parameter  weighting  matrix.  The  parameter 
trade-offs  that  result  from  singular  value  truncation  may  be  determined  by  an 
inspection  o'"  the  parameter  resolution  matrix, 

R  =  V  VT  ,  (16) 

which  reduces  to  the  identity  matrix  in  the  case  of  perfect  resolution.  If 


V  -  J*  V  V  V  V  W  V  •j<v> 


'J<  '/v  vvw  v?_v.v.v.v  v  v.-r.-7-,-.v; -’.v.  v.  j-. 


during  test  inversions,  the  parameter  resolution  is  deemed  unacceptable, 
adjustments  may  be  made  in  the  parameter  weights  or  the  singular  value 
cut-off,  and  the  inversion  may  be  repeated.  The  importance  of  particular 
observations  or  residuals  in  the  inversion  mny  be  determined  from  the 
information  density  matrix, 

S  =  U  UT  (]7) 


or  its  diagonal,  the  data  importance  vector.  This  may  indicate,  for  example, 
that  the  parameter  changes  are  most  sensitive  to  the  absolute  amplitude 
residual  at  nearer  receivers,  and  to  the  absolute  time  residual  at  more 
distant  receivers.  The  data  importance  vector  from  test  inversions  may 
suggest  modifications  to  the  observation  weights. 

Shaw  and  Orcutt  (1985)  use  a  different  approach  to  solve  the  linearized 
inverse  (which  they  attribute  to  Parker,  unpublished  manuscript,  1985)* 

Rather  than  inverting  the  equation 

Ac  =  A  Ap  (18) 

they  add  A  pg  (where  pg  contains  the  starting  model)  to  both  sides  and  invert 

(Ac  +  A  pQ)  =  A  (pQ  +  Ap)  =  A  pj  (19) 

to  solve  directly  for  the  new  model  parameters,  p-^.  The  benefit  of  (id), 
which  they  term  the  "Jumping"  method,  is  that  rather  than  minimizing  the 


parameter  changes  and  obtaining  a  model  dependent  on  the  starting  model,  one 


minimizes  (by  some  appropriate  norm)  the  resulting  model  itself.  An 
additional  benefit  is  that  while  a  priori  constraints  on  the  model  parameters 
(after  Jackson,  1979)  may  be  desired,  it  is  often  difficult  to  transform  these 
into  constraints  on  the  parameter  changes.  Although  we  have  not  yet 
implemented  or  tested  the  "jumping"  method,  it  would  seem  to  be  applicable  to 
the  near-field  inversion. 

TESTS  OF  THE  WAVEFORM  INVERSION  PROCEDURE 

Test  D  1  -  Ha If spa c e  Velocity 

As  a  simple  initial  tent  of  the  inversion  for  velocity  structure,  we 

generated  a  synthetic  data  set  for  an  explosion  source  in  a  half space.  The 

compressional  velocity  in  the  halfspace  was  6.0  km/s,  the  shear  velocity  was 

3.1*1*  km/s  (a  Poisson  solid)  and  the  density  was  3.0  g/cm^.  The  source  was 

represented  by  a  Helmberger-Hadley  (19fil)  RDP  with  K=7.0  s-1,  B=1.0  and 
10  7 

S’  *1.0x10  cm  ,  and  was  located  at  a  depth  of  1.0  km.  Vertical  component 
velocity  synthetics  were  computed  for  ranges  of  2,  U ,  6,  8,  10  and  12  km. 

Since  the  object  of  the  test  inversion  is  the  halfspace  velocity,  we  constrain 
the  structure  to  be  a  single  layer  with  a  zero  velocity  gradient. 

In  such  a  simple  case,  it  is  instructive  to  investigate  the  residuals  as 
a  function  of  the  parameter  values.  In  Figure  2  are  the  waveform  residuals, 
ei,  absolute  amplitude  residuals,  a^,  and  travel  time  residuals,  t^,  plotted 
as  functions  of  receiver  distance  and  halfspace  velocity.  Velocities  from  c . 0 
to  7.0  km/s  are  considered  In  increments  of  0.25  km/s.  Since  the  body 
waveforms  do  not  change  with  variations  in  the  halfspace  velocity,  the 
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e  2.  Waveform,  amplitude  and  travel-time  residuals  for  the  half spare 
model  of  Test  #  1  as  functions  of  receiver  distance  and  velocity.  The 
correct  velocity  is  6.0  km/s. 


variations  in  waveform  residuals  are  due  to  the  truncation  of  the  synthetic 
waveform  as  it  is  shifted  with  respect  to  the  time  window.  The  waveform 
residual  is  zero  only  when  there  is  a  perfect  match:  at  the  correct  halfspace 
velocity  of  6.0  kra/s.  Similarly,  the  absolute  amplitude  residual  does  not 
represent  the  difference  in  peak  amplitudes  between  observed  and  synthetic 
waveforms,  but  rather  the  difference  in  the  peak  values  of  their 
autocorrelations.  So,  systematic  amplitude  variations  due  to  velocity  changes 
are  superposed  with  truncation  effects  due  to  the  time  shifts  relative  to  the 
window.  Amplitude  residuals  may  take  on  positive  and  negative  values,  and  the 
sign  is  important  in  determining  the  partials.  Only  at  the  correct  velocity 
(6.0  km/s)  are  the  amplitude  residuals  at  all  of  the  receivers  zero.  The  most 
sensitive  measure  of  the  correct  halfspace  velocity  appears  to  be  the  travel 
time  residual.  This  residual  is  linearly  decreasing  at  all  receivers  for 
increasing  halfspace  velocity  and,  of  course,  the  most  distant  receivers  are 
most  sensitive  to  variations  in  the  velocity.  The  travel  time  residual  is 
zero  at  all  receivers  for  the  correct  velocity. 

For  simple  problems  with  only  a  couple  of  free  parameters,  a 
grid-searching  approach  in  which  the  minimum  of  the  residuals  is  found  by 
stepping  through  the  parameter  space  would  be  quite  reasonable,  particularly 
in  the  presence  of  local  minima.  This  is  essentially  the  procedure  used  to 
generate  Figure  2.  The  grid-searching  approach,  however,  becomes  quite 
unwieldy  when  the  number  of  free  parameters  exceeds  three  or  four,  and  does 
not  readily  supply  information  on  the  relative  resolution  of  the  parameters  or 
the  variance  of  the  solution.  In  anticipation  of  more  complex  problems, 
therefore,  we  have  performed  the  generalized  inverse  for  the  test  case  of 
determining  the  halfspace  velocity  from  a  starting  model  with  a  velocity  of 
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5.0  km/s.  If  we  were  to  invert  only  the  travel-time  residuals.  Figure  2 
indicates  that  the  inversion  would  be  linear.  However,  since  the  partials  are 
computed  by  a  linear  difference  of  the  residuals,  the  curvature  in  the 
waveform  and  amplitude  residuals  requires  an  iterative  scheme.  The  partials 
are  computed  for  perturbations  of  +10$  about  the  velocity  of  the  starting 
model.  Accounting  for  the  dimensionality  of  the  residuals,  and  utilizing  the 
observation  that  travel  time  residuals  appear  to  be  most  sensitive  to 
halfspace  velocity,  we  define  the  covariance  matrix  of  the  observations  to  be 
diag(l.O,  1.0,  0.1 )  for  the  waveform,  amplitude  and  time  residuals  at  each 
receiver,  respectively.  Since  there  is  only  one  parameter,  the  covariance  of 
the  velocity  is  set  to  1.0.  Five  iterations  were  allowed,  with  no  constraint 
or  damping  placed  on  the  parameter  changes.  The  results  are  listed  in  Table 
1,  and  the  waveforms  for  each  iteration  are  plotted  in  Figure  3. 

After  only  two  iterations,  the  velocity  determined  is  within  0.5$  of  the 
correct  value,  and  convergence  is  complete  after  four  iterations.  The  quality 
of  fit  may  be  determined  by  inspecting  the  waveforms  in  Figure  3,  or  by  two 
objective  measures  listed  in  Table  1.  The  RMS  fit  is  defined  as  RMS  = 

/r  /  3*ra,  where  r  is  the  objective  function  (11)  which  is  the  sum  of  the 
squared  residuals  and  3*ra  is  the  total  number  of  residuals.  The  least-squares 
error  is  LSE  =  |  A  Ap  -  Ac|^,  and  is  a  measure  of  the  linearity  of  the 
inversion  at  each  iteration.  In  initial  inversions,  the  time  windows  were 
chosen  to  include  the  halfspace  Rayleigh  wave  at  each  receiver.  This, 
however,  causes  a  strong  local  minimum  to  be  formed  at  low  velocities  when  the 
"synthetic"  P  wave  correlates  with  the  "observed"  Rayleigh  wave.  Although  the 
generalized  ray  synthetics  are  exact  for  the  halfspace  Rayleigh  wave,  they  are 
generally  inadequate  for  modeling  surface  waves  for  more  complicated 
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Table  1  -  Test  #  1  Results  -  Halfspace  Velocity 


Iteration 

a  (km/s) 

RMS  fit 

Least-square 

Start 

5.000 

0.277 

1 

5.629 

0.137 

1.978 

2 

5.969 

0.130 

0.280 

3 

5.987 

0.056 

0.308 

k 

5-999 

0.001 

0.056 

5 

5.999 

0.012 

0.0002 

Data 

6.000 

structures.  For  both  of  these  reasons,  the  inversion  time  windows  were 
shortened  to  exclude  the  Rayleigh  wave  where  it  is  separable  from  the  P  wave. 


Test  2  -  A  Single  Gradient 

For  a  somewhat  more  realistic  test,  and  in  order  to  test  the 
discretization  and  ray  summing  procedures,  a  synthetic  data  set  was  generated 
for  a  single  velocity  gradient.  The  corapressional  velocity  at  the  surface  was 
set  to  3-5  km/s  and  the  velocity  gradient  was  0.5  km/s/km.  The  shear  velocity 
was  constrained  to  that  of  a  Poisson  solid  (l//3  times  the  compressional 
velocity)  and  the  density  was  constrained  (for  simplicity  in  coding)  to  1/2 
the  compressional  velocity.  The  gradient  was  discretized  using  a  predominant 
period  of  0.25  s,  which  results  in  10  layers  in  the  top  15  km  with  layer 
thickness  increasing  with  depth.  All  upgoing  and  diving  P  waves  and  the 
surface-reflected  pP  waves  were  considered,  but  culling  late  arriving  rays 
resulted  in  18  rays  being  summed  for  the  final  response.  Once  again. 
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synthetic  velocity  seismograms  were  computed  for  ranges  of  2,  4,  6,  fi,  10  and 
12  km. 

Figures  4  and  5  illustrate  the  dependence  of  the  waveform,  absolute 
amplitude  and  travel  time  residuals  at  each  receiver  as  functions  of  the 
velocity  gradient  and  the  velocity  at  the  surface,  respectively.  As  shown  in 
Figure  4,  since  the  waveforms  at  the  nearer  receivers  are  dominated  by  the 
upgoing  P  wave,  which  is  relatively  insensitive  to  the  velocity  gradient, 
there  is  very  little  variation  in  their  residuals  over  the  ranges  of  gradient 
considered  (0.3  to  0.7  km/s/km).  On  the  other  hand,  the  waveforms  at  the 
farther  receivers  are  dominated  by  diving  turning  rays,  and  their  residuals 
(particularly  travel  time)  resolve  changes  in  the  gradient  quite  well. 

Overall,  it  appears  that  velocity  gradient  in  this  simple  example  will  be 
resolved  by  travel  time  residuals  at  farther  receivers  and  amplitude  residuals 
at  all  receivers.  The  residuals  due  to  variations  in  surface  velocity  (Figure 
5)  are  very  similar  to  those  for  the  halfspace  test  (Figure  2).  The  waveform 
residuals  have  a  broad  minimum,  the  travel  time  residuals  are  a  robust  measure 
of  the  correct  surface  velocity,  and  the  amplitude  residuals  include 
variations  due  to  truncation  by  the  time  window.  In  comparing  Figures  4  and 
5,  it  appears  that  the  residuals  are  more  sensitive  to  variations  in  surface 
velocity  than  velocity  gradient.  We  will  have  to  take  this  into  account  in 
assigning  parameter  weights  in  the  inversion. 

Before  performing  an  inversion  for  surface  velocity  and  gradient,  it  is 
instructive  to  invert  for  gradient  only,  constraining  the  surface  velocity  to 
its  correct  value.  In  this  case,  we  define  the  covariance  matrix  of  the 
observations  according  to  the  dimensionality  of  the  residuals,  but  in 
addition,  we  wish  to  preferentially  minimize  the  waveform  residuals  since 
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e  U.  Waveform,  amplitude  and  travel-time  residuals  for  the  gradien 
model  of  Tests  #  2  and  3  as  functions  of  receiver  distance  and  velo 
gradient.  The  correct  gradient  is  0.5  km/s/km. 
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Figure  5*  Waveform,  amplitude  and  travel-time  residual 
model  of  Tests  P  2  and  3  as  functions  of  receiver 
velocity  at  the  surface.  The  correct  velocity  is 
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changes  in  gradient  appear  to  have  significant  effect  on  the  details  of  the 
waveform.  We  define  the  matrix  to  be  diag(0.01,  1.0,  O.l)  for  waveform, 
amplitude  and  time  residuals  at  each  receiver,  respectively.  Partial 
derivatives  are  computed  by  perturbations  in  the  gradient  of  ^10^,  beginning 
with  a  gradient  of  0.6  km/s/km.  The  correct  value  is  0.5  km/s/km.  Five 
iterations  were  allowed  with  no  damping  of  parameter  changes.  The  results  of 
the  inversion  are  listed  in  Table  2  and  the  waveforms  for  each  iteration  are 
plotted  in  Figure  6. 

Convergence  to  the  solution  is  much  more  difficult  in  this  case.  An 
inspection  of  the  data  importance  vectors  (Figure  7)  indicates  that,  the  first 
iteration  is  primarily  controlled  by  the  waveform  and  travel  time  residuals  at 
the  12  km  receiver,  with  some  influence  by  the  waveform  and  amplitude 
residuals  at  6  km.  After  this  initial  step  toward  the  correct  solution,  the 
second  and  third  iterations  are  controlled  by  the  amplitude  residuals  at  6  and 
10  km.  As  may  be  seen  in  Figure  U,  the  amplitude  residuals  at  these  receivers 
vary  more  rapidly  for  gradients  slightly  less  than  0.5  km/s/km  than  for  those 
slightly  greater  than  0.5  km/s/km.  In  the  second  iteration,  the  partial 
derivatives  are  computed  based  on  the  difference  in  the  residuals  at  0.573  and 
0.469  km/s/km.  Even  though  the  latter  is  closer  to  the  correct  solution  (0.5 
km/s/km),  its  amplitude  residual  is  larger  in  absolute  value  than  the  former, 
so  the  partial  derivative  of  the  amplitude  residuals  causes  the  solution  to 
move  away  from  the  correct  gradient.  By  the  fourth  iteration,  the  solution 
has  noved  far  enough  away  that  *he  lower  bracket  value  used  in  computing  the 
partials  is  approximately  0.5  kn/s/km,  and  the  partial  is  unaffected  by  the 
changes  in  the  residuals  for  smaller  gradients.  Second  derivative 
information,  available  by  considering  the  residuals  of  the  starting  model  as 


Table  2  -  Test  #  2  Results  -  A  Single  Gradient 


Iteration 

Vet  (km/s/km) 

RMS  fit 

Least-square 

Start 

0.600 

0.102 

1 

0.522 

0.052 

1.015 

2 

0.529 

0.107 

0.125 

3 

0.560 

0.083 

2.12R 

1* 

O.U97 

0.01*7 

0.080 

5 

0.1*71 

0.117 

0.088 

Data 

0.500 

well  as  the  bracket  values,  could  minimize  this  instability.  However,  in  an 
effort  to  reduce  computational  effort  in  inversions  for  larger  numbers  of 
parameters,  partial  derivatives  will  be  computed  from  a  one-sided  perturbation 
from  the  starting  model.  In  addition  to  nearly  halving  the  required 
computations,  this  also  decreases  the  range  of  parameter  space  over  which  the 
linear  approximation  for  partials  is  made.  The  adequacy  of  the  linear 
approximation  is  quantified  by  the  least-squares  error.  Often  the 
least-squares  error  will  actually  increase  as  the  correct  solution  is 
approached  due  to  the  curvature  of  the  residuals  about  that  value. 


Test  H  3  -  Velocity  and  Gradient 

We  have  performed  several  test  inversions  for  simultaneously  determining 
velocity  gradient  and  surface  velocity.  There  is  an  obvious  trade-off  between 
these  parameters,  indicating  the  need  to  truncate  small  singular  values.  In 
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addition,  however,  the  parameters  should  be  weighted  beyond  that  required  by 
the  dimensionality  of  the  parameters.  For  the  synthetic  data  set  of  the 
previous  test.  Figures  U  and  5  illustrated  that  the  inversion  appears  to  be 
more  sensitive  to  variations  in  surface  velocity  than  gradient.  This  means 
that  a  given  residual  value  may  be  due  either  to  a  reasonable  variation  in 
surface  velocity  or  a  relatively  larger  variation  in  gradient.  Since  the 
inversion  is  non-linear,  an  error  in  surface  velocity  will  be  only  partially 
corrected,  while  the  gradient  varies  substantially  to  compensate.  Inversions 
in  which  the  parameters  are  equally  weighted,  or  in  which  gradient  is 
preferentially  weighted  will  diverge.  For  this  test,  we  have  therefore 
defined  the  covariance  matrix  of  the  parameters  such  that  surface  velocity  is 
preferentially  weighted  as  diag(l.O,  0.1 )  for  surface  velocity  and  gradient, 
respectively.  The  covariance  matrix  of  the  observations  is  diagfO.l,  1.0, 

0.1)  for  waveform,  amplitude  and  travel  time  residuals,  respectively.  A 
cut-off  was  established  so  that  one  of  the  singular  values  could  be  truncated 
to  avoid  large  variations  due  to  poor  resolution.  Perturbations  of  +10^  were 
used  to  generate  the  partials.  Ten  iterations  were  allowed  from  a  starting 
nnb'l  with  sirfaee  velocity  of  l?.1^  km/s  and  a  gradient  of  0.U  km./s/km.  The 
correct  solution  h  id  a  v*» »-ity  of  3.5  km/s  and  a  gradient  of  0.5  km/s/km. 

"’•ip  r*«.:  .  >  r  *  v*  inversion  are  listed  in  Table  3,  the  waveforms  from 

:v':  in  vigure  °,  and  the  resolution  matrix  for  each 
it  era*.;  n  .  ■  k  :g  ;r-  K  In  the  first  iteration,  both  parameters  are 

perfer*  y  r  >-  ..  .1  1  *  n~  'iv.,'"  i  s  toward  the  correct  solution  for  both.  In 

the  se  ■  .  i  -i  :  •  <w .  •  .•  one. ,  one  singular  value  has  been  truncated. 

mhp  weigv  .  •  h-i  ;  ,.e<  m,h*  the  remaining  eigenvector  is  oriented  in 

th"  pa r  i.m> •  •  ^  *  ly  ii:  •  n<-  -direct  i  >» ,  ,,«•  surface  velocity,  but  with  a 


Table  3  -  Test  #  3  Results  -  Velocity  and  Gradient 


Iteration 

ft  (km/s) 

R* 

Vix  (km/s/km) 

R* 

RMS  fit 

L3E 

Start 

2.500 

0.1*00 

0.598 

1 

2.751 

1.000 

0.1*95 

1.000 

0.311 

25.9 

2 

2.991* 

0.967 

0.509 

0.033 

0.238 

10.1 

3 

3.331 

0.921 

O.5UO 

0.079 

0.156 

1.01 

1* 

3.»*53 

0.936 

0.550 

0.061* 

0.139 

0.1*09 

5 

3.1*58 

0.81*2 

0.551 

0.158 

0.091* 

0.369 

6 

3.1*37 

0.91*7 

0.51*9 

0.053 

0.11*1* 

0.181* 

7 

3.1*39 

0.91*7 

0.51*9 

0.053 

0.151 

0.1*06 

8 

3.1*32 

0.917 

0.51*9 

0.083 

0.11*0 

0.1*70 

9 

3.1*56 

0.967 

0.550 

0.033 

0.129 

0.363 

10 

3.1*66 

0.881 

0.551 

0.119 

0.078 

0.311 

Data 

3-500 

0.500 

value  of  the  diagonal  of  the  resolution  natrix  for  this  parameter. 


small  component  in  the  gradient.  This  means  that  the  second  and  following 
iterations  consist  mainly  of  modifications  to  the  surface  velocity,  with  only 
small  (and  perhaps  incorrect)  modifications  to  the  gradient.  After  about  five 
iterations,  the  model  is  slightly  low  in  surface  velocity,  which  is 
counterbalanced  by  a  slightly  high  gradient.  The  waveforms  from  this 
combination  are  close  enough  to  those  of  the  correct  solution,  that  parameter 
perturbations  are  relatively  small  in  later  iterations.  Figure  8  shows  that 
the  primary  effect  of  these  later  perturbations  is  to  improve  the  absolute 
amplitude,  since  travel  time  and  waveshape  are  already  fairly  well  fit.  This 
indicates  that  it  mny  be  desirable  to  allow  the  covariance  matrix  of  the 
observations  to  change  for  later  iterations,  to  allow  the  amplitude  residual 
more  influence  on  the  solution. 


I-  3? 


I"  34 


x 

P  4J 
to 

>«  4-5 
•H  gJ 

v-  E  * 

CO 

U  s  C 
.c  v  o 

P  *H  P 
P  P 

c  c  cfl 

H  d'  V- 

T3  a; 

*H  .p 
•  »«H 

m  a? 

.c  bO 
=»t  p  c 

•H 

P  CO  > 

W  *rH  C 
0)  fP 
^  X  P 

•H  C 

dn  ^  P 
O  P 

c  6,5 

o 

•p  a) 

p  ,c  • 

Cl!  P  '— > 

p  .c 

1’  >C  P  • 
P  C  M  « 
•h  nj  c  aj 
<D  >■ 

r  'Cr'rl 

O  0)  o 
Ct)  >  P  CO 

l)  H  .rl  H 
CEP 

p  to  3 

O  B  P 

ft  p  a>  a> 
>  p 
to  >>  a)  P 
D  ftr  ,c  a; 
op  ,o 
•p  o  s, 
pops, 

p  tp  -h  +? 

aj  p  o  -p 
E  o  O  V 
ftrP  O 
C  r-l 

o  o  >•  o 

-p  p  >■ 
p  ((« 

3  C  £ 

-p  to  (C  P 
O  P  *H 
m  t  p  It 
o  p  c 

P  fl)  o  * 
E  -P  ft 
P  et!  T3  ft 
0  p  nj  C 
P  gi  p 
|  ft  00  ti 

ttf  O  J3  ’S 
P  ,C  P  P 
rf  P  O  P 
ft  ,£> 

-  CO 
O  C  P  P 

•coon 

E-*  -p  C  p 
P  4) 
rf  to  fc 

•  p  y  ec$ 
o\  o  o  p 
P  -P  ctf 
O  *P  p  ft 
p 


J«  '-r  ’*4  *,* 


y-  y  t*'  ysy*  yyr  yv 


Test  Source  and  Structure 

Finally,  in  order  to  test  the  simultaneous  inversion  for  source  and 
structure  parameters,  we  generated  a  synthetic  data  set  using  a 
Helraberger-Hadley  (’981)  RDP  and  a  two-layer  gradient  structure.  The  velocity 
structure  was  based  on  the  model  Hartzell,  et  al.  (1983)  obtained  for  Pahute 
Mesa,  and  consists  of  a  surface  velocity  of  2.5  km/s  and  gradient  of  O.83 
km/s/km  in  the  top  layer.  The  second  layer  is  separated  from  the  first  by  a 
second-order  discontinuity  at  a  depth  of  3.0  km  and  has  a  gradient  of  0.20 
km/s/km  below  that.  For  the  source  function,  we  set  ¥  =  1.0  x  10^°  cm^, 

OO 

k=7.0  s-"'-  and  B=l.  Vertical  component  velocity  synthetics  were  computed  for  a 
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source  at  1.0  km  depth  and  for  ranges  from  2  to  12  km,  later  extended  to  20 
km. 

Figure  10  shows  the  residuals  for  the  six  free  parameters  of  a  two-layer 
model  with  a  second-order  discontinuity  for  ranges  of  2-20  km.  In  the  figure, 
all  of  the  residuals  have  been  plotted  on  a  scale  with  a  maximum  amplitude  of 
1.0  so  that  the  relative  resolution  of  the  parameters  may  be  inferred  at  a 
glance.  4^  affects  only  the  amplitude  residual,  in  a  linear  fashion  and,  of 
course,  independent  of  range,  k  is  also  independent  of  range  and  affects 
predominantly  the  amplitude  residual,  but  values  well  away  from  the  correct 
one  (7.0  s“M  also  affect  the  waveform  and  travel  time  residuals.  This  is 
because  a  change  in  rise  time  of  the  source  affects  the  interference  phenomena 
that  determine  the  waveform  and  time  of  maximum  correlation.  Obviously  4'  and 
k  will  trade  off,  but  k  also  includes  trade-offs  with  other  parameters.  By 
far  the  most  important  parameter  for  the  waveform  and  travel  time  residuals  is 
the  velocity  at  the  top  o1’  the  first  layer.  As  in  the  halfspace  velocity 
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Figure  10.  Waveform,  amplitude  and  travel-time  residuals  for  the  six 

inversion  parameters  of  Test  0  I*  as  functions  or  receiver  distance.  All 
the  plots  have  boon  scaled  to  a  maximum  residual  value  of  1.0,  so  tnat 
the  relative  sensitivity  of  the  residuals  to  various  parameters  may  be 
Inferred  at  a  glance.  The  correct  value  for  each  parameter  is  located  at 

the  center  of  each  plot.  ^ 
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test,  the  travel  time  residual  is  linearly  related  to  top  velocity  and  the 
effect  increases  with  range  until  portions  of  the  waveform  fall  outside  of  the 
time  window.  As  in  Test  #  2,  the  upper  gradient  is  hardly  resolved  at  all  at 
the  nearest  receivers,  and  the  residual  minima  are  quite  broad  even  at  the 
farthest  ranges.  Due  to  the  depth  of  the  second  layer  in  our  synthetic  model, 
the  lower  gradient  is  unresolved  out  to  10-12  km  and  the  residuals  are 
relatively  quite  small  out  to  20  km.  Similarly,  the  interface  depth  parameter 
is  resolved  only  when  it  approaches  the  source  depth  (l.O  km).  Of  course,  if 
the  top  layer  parameters  are  near  their  correct  values,  the  residuals  will  be 
relatively  more  sensitive  to  variations  in  lower  layer  parameters. 

A  reasonable  procedure  for  this  and  following  inversions  is  to  start  with 
a  simple  parameterization  to  obtain  preliminary  estimates,  then  to  increase 
the  complexity  of  the  model  and  refine  the  parameter  estimates.  This  is 
particularly  advisable  when  certain  parameters  are  far  better  resolved  than 
others.  For  example,  we  could  solve  initially  for  4'  and  k  (holding  B 
constant)  and  the  parameters  of  a  single  gradient.  Using  these  as  starting 
values,  we  could  then  assume  a  second-order  discontinuity  and  attempt  to 
resolve  the  gradient  and  depth  of  the  second  layer.  In  addition  to  source  and 
upper  layer  parameters,  we  are  then  inverting  for  a  total  of  six  parameters. 
Following  this,  we  could  free  the  second  layer  to  have  a  first-order 
discontinuity  (include  the  velocity  at  the  top  of  the  second  layer),  solving 
for  seven  parameters.  As  the  estimates  of  the  better  resolved  parameters 
approach  their  correct  values,  we  could  include  a  third  or  even  a  fourth 
layer,  consider  anelastic  attenuation,  or  allow  B  to  vary. 

Table  U  and  Figure  11  include  the  results  of  the  first  two  steps  of  this 
inversion  procedure  applied  to  the  synthetic  data  set.  In  the  first  five 
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Table  *4  -  Test  tt  4  Results  -  Source  and  Structure 


Iteration 

H’oo 

( xlO^cnP ) 

k 

(s'1) 

VOL 

(km/s/km) 

(km}s ) 

Vou 

(km/s/km) 

7 

(km) 

RMS  fit 

Start 

2.00 

9.00 

0.50 

2.00 

1.098 

1 

1.33 

6.70 

0.63 

2.50 

- 

- 

0.298 

2 

1.33 

6.71 

0.78 

2.51* 

- 

- 

0.231 

3 

1.33 

6.7? 

0.83 

2.50 

- 

- 

0.253 

4 

1.33 

6.70 

0.87 

2.U7 

- 

- 

0.245 

5 

1.3? 

t>.6H 

0.95 

2.30 

- 

- 

0.228 

[Startl 

11.27] 

16.371 

1  10.89] 

I2.U5] 

10.89] 

12.00] 

[0.2991 

6 

0.72 

6,20 

0.89 

2.46 

0.58 

2.?9 

0.208 

7 

0.814 

6.2? 

0.89 

2.44 

0.49 

2.56 

0.200 

8 

0.9)4 

6.25 

0.89 

2.49 

0.44 

2.57 

0.174 

9 

0.91 

0.2  >4 

0.89 

2.44 

0.43 

2.57 

0.183 

10 

0.91 

6.24 

0.89 

2.45 

0.39 

2.68 

0.165 

Data 

1.00 

7.00 

0.83 

2.60 

0.20 

3.00 

Jf. 

The  starting  model  for  iterations  6-10  was  inadvertantly  changed 
from  the  model  of  iteration  5* 


iterations,  we  solve  for  the  source  and  upper  layer  parameters  using  synthetic 
data  at  ranges  of  2-12  km.  The  parameter  weighting  matrix  was  defined  as 
diag(0.01,  0.10,  0.01,  1.00)  for  k,  gradient,  and  top  velocity, 
respectively.  The  covariance  matrix  of  the  residuals  was  diag(0.10,  0.10, 
0.001)  for  waveform,  amplitude  and  travel  time  residuals  at  receivers  from  2-6 
km,  respectively.  In  order  to  improve  resolution  of  the  gradient,  the 
waveform  covariance  of  receivers  from  8-1?  km  was  decreased  so  that  the 
covariance  matrix  became  dlag(0.01,  1.00,  0.001).  Figure  12  presents  the 
parameter  resolution  matrices  for  each  iteration,  and  shows  that  the  upper 
layer  parameters  are  perfectly  resolved  through  the  first  five  iterations.  ' , 
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:  parameter  resolution  matrices  for  each  iteration  of  Test  *  U.  From  top  to  bo* 
or respond  to  the  source  parameters,  Tm  and  k,  and  the  structure  parameters,  upp 
top  velocity,  lower  gradient  and  interface  depth.  In  iterations  1-5,  the  lower 
parameters  are  not  included  in  the  inversion. 


and  k  trade-off  equally  in  the  first  iteration  and  are  unresolved  in  later 
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iterations.  For  the  most  part,  the  data  importance  vectors  (Figure  13) 
reflect  the  residual  covariances.  In  the  first  iteration,  the  solution  is 
controlled  by  the  travel  times  and  amplitudes  of  the  nearest  receivers,  while 
in  iterations  2-5,  waveforms  and  travel  times  to  the  farther  receivers  becomes 
important . 

Several  attempts  were  made  to  incorporate  the  lower  layer  parameters  into 
the  inversion  using  the  results  of  iteration  5  as  a  starting  model.  However 
using  only  receivers  from  2-12  km,  the  lower  gradient  and  depth  were  poorly 
resolved  and  the  resulting  inversions  were  unstable.  At  this  point,  we 
extended  the  synthetic  data  set  by  computing  waveforms  for  receivers  at  lU-20 
kin.  The  covariance  matrix  of  the  residuals  at  all  receivers  was  set  as 
diag(0.10,  0.10,  0.001)  and  the  parameter  weighting  matrix  was  diag(0.10, 

0.10,  0.001,  0.01,  1.0,  1.0)  for  k,  upper  gradient,  top  velocity,  lower 
gradient  and  interface  depth,  respectively.  The  discontinuity  was  constrained 
to  second  order. 

As  shown  in  Table  U,  rather  than  using  the  results  of  Iteration  5  as  the 
starting  model  for  iterations  6-10,  the  results  of  an  earlier  inversion  were 
inadvertantly  used.  The  model  parameters  are  quite  similar,  but  the  weighting 
scheme  of  the  earlier  run  resulted  in  poorer  convergence  (see  RMS  fit). 
Ideally,  the  starting  model  should  be  irrelevant  to  the  final  solution,  but 
due  to  nonlinearity  and  the  parameter  weighting  designed  to  stabilize  the 
inversion,  this  is  not  entirely  the  case.  The  parameter  resolution  matrices 
(Figure  12)  show  that  while  and  the  lower  gradient  and  depth  are  well 
resolved,  the  upper  gradient  and  k  are  poorly  resolved.  The  top  velocity  is 
resolved,  but  trades  off  with  the  upper  gradient.  The  results  for  iterations 
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6-10  (Table  U),  therefore,  indicate  that  while  'l,uo,  the  lower  gradient  and  the 
depth  approach  their  correct  values  (shown  at  the  bottom  of  the  table),  the 
parameters  k,  upper  gradient  and  top  velocity  depart  slightly,  if  at  all,  from 
the  starting  model.  In  the  case  of  the  top  velocity,  this  is  because  the 
starting  model  for  iteration  6  is  quite  near  the  correct  value  (the  correct 
value  was  obtained  in  the  first  iteration  of  the  inversion).  The  data 
importance  vectors  (Figure  13)  show  that  in  iterations  6-10,  the  travel  time 
and  amplitude  of  the  farther  receivers  are  most  important  in  determining  the 
solution,  although  the  travel  time  at  near  receivers  is  important  as  a 
constraint.  Waveform  residuals  contribute  little  to  the  solution.  This  may 
also  be  seen  in  the  plots  of  the  waveforms  (Figure  11 ),  in  which  the  character 
of  the  data  has  been  well  matched  by  the  fifth  iteration,  but  the  travel  time 
and  amplitude  of  farther  stations  needs  refinement. 

The  results  of  this  test  inversion  may  be  summarized  by  the  structure 
models  that  are  obtained  (Figure  lU).  In  this  plot  of  P-wave  velocity  vs. 
depth,  the  model  from  which  the  synthetic  starting  model  is  shown  as  the  broad 
dashed  line  labeled  "START".  After  five  iterations,  the  top  velocity  and 
upper  gradient  have  been  matched  almost  exactly  (lower  structure  was  not 
included  in  these  iterations).  The  RDP  parameters  reached  their  best 
estimates  for  this  parameterization  within  the  first  iteration  and  remained 
unchanged  to  iteration  five.  The  starting  model  for  iterations  6-10  was 
similar  to  the  iteration  5  solution.  The  lower  gradient  was  assumed  to  be  the 
same  as  the  upper  gradient  and  the  depth  to  the  top  of  the  second  layer  was 
assumed  as  2  km.  The  details  of  the  structure  between  2  and  3  km  could  not  be 
resolved  with  receivers  only  out  to  12  km,  so  the  data  set  was  extended  to 
ranges  up  to  20  km.  In  iterations  6-10,  the  upper  layer  parameters  remained 


Test  $  4  Inversion  Structure  Models 


P  VELOCITY  (km/s) 


2.0  3.0  4.0  5.0  0.0  7.0 


Figure  lU.  Velocity  structure  models  from  Test  §  *4.  The  synthetic  "data"  set 
was  computed  for  the  two-gradient  model  plotted  as  a  solid  line.  The 
starting  model  consisted  of  a  single  gradient  shown  as  the  long  dashed 
line.  In  the  first  five  iterations,  a  single  gradient  was  considered, 
and  the  resulting  model  is  denoted  "It.  5"»  It  matches  the  upper 
gradient  of  the  "data"  model  quite  well.  Finally,  a  two-gradient  model 
was  obtained  after  iteration  10  (short  dashed  line)  which  approaches  the 
lower  "data"  gradient. 
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essentially  unchanged,  the  depth  increased  toward  the  correct  value  of  3.0  km 
and  the  lower  gradient  decreased  toward  the  correct  value.  The  source 
parameters  adjusted  mostly  in  the  sixth  and  seventh  iterations,  then  remained 
stable  through  the  remaining  iterations. 
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We  begin  the  inversions  of  actual  data  with  the  Pahute  Mesa  event,  BOXCAR 
(U/26/68,  1300  kT  announced).  Hartzell,  et  al.  (1983)  developed  their  Pahute 
Mesa  structure  model  primarily  by  modeling  the  waveforms  of  receivers  S-2l*  and 
S— 3U  from  BOXCAR,  although  travel  times  from  several  other  events  were  quite 
well  predicted.  We  will  use  an  approximation  of  the  Hartzell,  et  al.  (1983) 
model  (the  "data"  model  from  test  #  U)  as  our  starting  model  and  invert 
directly  for  a  two-layer  model  with  a  second  order  discontinuity.  The 
stations  used  in  the  inversion  are  listed  in  Table  5,  and  include  a  broad 
distribution  of  ranges.  The  closest  receivers  (S-12  at  3.81  km  and  S-l6  at 
U.87  km)  are  only  marginally  outside  the  spall  zone  and  appear  to  include 
arrivals  attributable  to  spall  opening  or  slap-down.  Since  these  effects  will 
not  be  included  in  our  modeling,  we  must  increase  the  covariance  of  the 
waveform  residuals,  effectively  downweighting  the  importance  of  fitting  the 
waveform,  at  these  receivers.  The  covariance  matrix  is  diag(l.00,  0.01, 

0.001)  for  waveform,  amplitude  and  travel  time,  respectively.  Similarly,  the 
farthest  station  (S-7^  at  22, U7  km)  is  located  outside  the  Silent  Canyon 
Caldera  of  Pahute  Mesa,  and  we  could  expect  the  plane-layered  structure 
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Table  5  -  Pahute  Mesa  Station  Parameters 


Station  Range  Window  Start  Window  End 
(km)  (sec)  (sec) 


BOXCAR 


INLET 


MAST 


SCOTCH 


S-12 

3.81 

S-l6 

U.B7 

S-2U 

7-27 

S-31* 

10.37 

S-71* 

22.  1*7 

S-5 

1.63 

S-6 

3.27 

S-7 

6.53 

S-5 

3.65 

S-6 

5-U7 

S-7 

7.30 

S-3A 

U.13 

S-k 

6.06 

1.25 

2.05 

1.55 

2.35 

1.85 

3.20 

2.70 

3.95 

1+.  50 

6.00 

o.uo 

1.25 

0.90 

1.75 

1.85 

2.70 

1.00 

1.90 

1.1*0 

2.35 

1.85 

2.80 

1.20 

2.00 

1.75 

2.30 

2.70 

** 

### 


Time  after  origin. 

End  time  for  iterations  1-5. 
End  time  for  Iterations  6-15. 


approximation  to  break  down  for  this  receiver.  We  have  therefore  increased 
the  covariance  of  the  waveform  and  travel-time  residuals  at  this  receiver. 

The  resulting  covariance  is  diag(0.10,  0.01,  0.10).  This  leaves  two  receivers 
(S-2U  at  7.27  km  and  S-3^  at  10.37  km)  with  reasonable  covariance  in  all  three 
residuals.  The  covariance  for  these  receivers  Is  diag(0.01,  0.01,  0.001 ). 
After  many  test  inversions,  the  parameter  weighting  scheme  that  allowed  a 
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stable  solution  was  found  to  be  diag(0.10,  0.10,  0.001,  0.01,  0.01,  0.10)  for 
^oo,  k,  upper  gradient,  top  velocity,  lower  gradient  and  interface  depth, 
respectively.  Five  iterations  were  allowed,  truncating  one  eigenvalue  in  each 
iteration. 

The  inversion  results  are  listed  in  Table  6  and  the  waveforms  are 
displayed  in  Figure  IS.  The  synthetic  waveforms  for  the  nearest  receivers  are 
simple  and  impulsive,  while  the  observed  waveforms  are  broad,  and  include 
significant  late  arrivals.  Since  the  travel-time  residual  is  defined  as  the 
peak  in  cross-correlation  with  the  time  window  denoted  by  arrows  in  Figure  15, 
the  minimum  of  the  residual  corresponds  to  the  alignment  of  peaks  and  not 
necessarily  the  alignment  of  f'rst  breaks.  This  error  due  to  unmodeled 
arrivals,  combined  with  delays  due  to  differing  receiver  structure  (considered 
noise)  accounts  for  the  0.17  sec  travel-time  error  at  the  7  km  station.  We 
expect,  and  have  allowed,  larger  travel-time  errors  at  the  farthest  receiver, 
although  the  final  synthetic  is  only  0.31  sec  late.  Although  the  numbers  to 
the  right  of  the  waveforms  in  Figure  15  are  the  maxima  of  the  trace  absolute 
values  and  not  the  first  peak  amplitudes,  we  see  that  the  final  synthetic 
first  peak  amplitudes  are  in  reasonable  agreement  with  the  observations  at  all 
but  the  nearest  receivers.  Details  of  the  waveforms  of  the  farther  three 
receivers  are  quite  well  modeled,  even  to  a  certain  extent  outside  the  time 
window  of  the  inversion.  The  first  peak  at  the  7  and  10  km  stations  includes 
the  constructive  interference  of  upgoing  and  diving  rays,  while  the  first 
trough  includes  the  interaction  of  the  surface  reflection  phases.  The 
waveform  at  the  farthest  station  changes  significantly  through  the  iterations 
in  an  attempt  to  fit  the  double  first  peak  and  the  gradually  rising  trough. 
Although  time  shifted,  the  synthetic  from  the  fifth  iteration  provides  a 
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Table  6  -  BOXCAR  Inversion  Results 


Iteration 

too 

(xlO^cm^) 

k 

(s"1) 

VOL 

(km/s/km) 

(km}s ) 

Votp 

(km/s/km 

Z2 
)  (km) 

RMS  fit 

Start 

12.00 

6.50 

0.83 

2.50 

0.20 

3.00 

0.2965 

1 

12.52 

(0.20) 

7.96 

(0.69) 

0.93 

(0.30) 

2.1+7 

(0.16) 

0.11 

(0.12) 

3.26 

(0.28) 

0.1919 

2 

12.1+6 

(0.1U) 

7.33 

(0.1+9) 

0.77 

(0.23) 

2.76 

(0.32) 

0.13 

(0.0M 

3.15 

(0.1+0) 

0.1999 

3 

12.56 

(0.15) 

7.51 

(0.1+3) 

0.85 

(o.uu) 

2.1+7 

(0.21) 

0.12 

(0.07) 

3.18 

(0.72) 

0.2313 

k 

12.78 

(0.20) 

8.06 

(0.51) 

0.88 

(0.11) 

2.1+3 

(0.16) 

0.15 

(0.01+ ) 

3.22 

(0.27) 

0.2615 

5 

12.63 
(0.1U ) 

7.62 

(0.1+2) 

0.88 

(0.2M 

2.U7 

(0.17) 

0.16 

(0.05) 

3.60 

(0.60) 

0.1801 

—  Numbers  in  parentheses  are  the  standard  deviations  of  the  parameter 
changes  that  result  in  the  model  parameters  for  that  iteration. 


reasonable  fit  to  these  features. 

Table  6  includes,  in  parentheses,  the  standard  deviations  of  the 
parameter  changes  for  each  iteration.  As  discussed  previously,  these  are  not 
the  standard  deviations  of  the  parameters  themselves,  and  could  only  be 
interpreted  so  if  the  starting  model  of  the  iteration  was,  in  some  sense, 
error  free.  This  error  measure  also  includes  the  effects  of  parameter  and 
residual  weighting  for  stabilizing  the  non-linear  inversion.  While  the 
standard  deviations  seem  to  reflect  mostly  the  sizes  of  the  parameter  changes 
in  a  given  iteration,  they  indicate  at  least  the  level  of  confidence  that  my 
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3.81  km 


65.043 

105.250 

89.668 

91.259 

103.830 

104.150 

87.677 


7.27  km 


UW 


43.469 

41.773 

34.578 


40.224 

37.112 

40.850 


22.47  km 


i 

m 


4.87  km 


10.37  km 


_J\ 


65.032 

61.710 

66.835 

67.052 

59.666 

76.903 

99.556 


14.406 

18.747 

17.282 


17.756 


18.281 

21.461 

21.931 


Figure  15.  Inversion  results  for  the  Pahute  Mesa  event,  BOXCAR.  Shown  are 
synthetic  seismograms  for  the  starting  model  and  each  iteration  (above) 
compared  with  the  observed  waveforms  (below).  The  peak  trace  amplitude 
is  shown  to  the  right  and  the  inversion  time  window  is  denoted  by  arrows 
The  time  scale  is  shown  on  the  bottom  right. 
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be  placed  on  a  particular  parameter  obtained  in  the  inversion.  The  parameter 
resolution  matrices  for  the  BOXCAR  inversions  are  plotted  in  Figure  l6.  The 
structure  parameters  are  perfectly  resolved  through  all  iterations.  4^  and  k 
trade  off,  with  k  being  somewhat  better  resolved.  The  data  importance  vectors 
(Figure  IT)  reflect  primarily  the  residual  weighting  scheme  discussed  above, 
but  indicate  that  waveforms  at  farther  stations,  travel  times  to  the  nearest 
stations  and  amplitudes  at  all  stations  influence  the  solution. 

The  final  structure  model  for  BOXCAR  is  not  significantly  different  than 
the  two-gradient  starting  model  based  on  the  model  of  Hartzell,  et  al.  (1983). 
Both  4'  and  k,  however,  are  significantly  larger  than  the  starting  model 
values,  which  were  based  on  the  values  obtained  by  Barker,  et  al.  (1985)  using 
the  Hartzell,  et  al.  (1983)  model.  In  addition,  the  RMS  fit  has  been 
significantly  improved  over  the  starting  model,  even  in  the  presence  of  noise 
and  unmodeled  arrivals.  An  inversion  of  the  BOXCAR  data  set,  in  which  the 
depth  of  the  water  table  was  fixed  as  a  layer  interface,  and  the  shallow 
structure  model  of  Leonard  and  Johnson  (1986)  was  used  as  the  starting  model, 
resulted  in  an  intermediate  layer  with  vanishing  thickness,  and  relatively 
poor  waveform  fits.  This  indicates  that  the  details  of  structure  above  the 
source  are  difficult  to  resolve  uniquely  from  waveforms  located  outside  the 
spall  zone. 
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BOXCAR  Resolution  Matrices 
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Figure  l6.  Parameter  resolution  matrices  for  each  iteration  of  the  BOXCAR 
inversion.  Plotted  from  top  to  bottom  are  the  rows  of  the  matrices 
corresponding  to  the  source  parameters,  'Vm  and  k,  and  the  structure 
parameters,  upper  gradient,  top  velocity,  lower  gradient  and  interface 
depth.  For  perfect  resolution,  the  resolution  matrix  would  be  the 
identity  matrix. 
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BOXCAR  Data  Importance  Vectors 
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Figure  IT.  Data  importance  vectors  for  each  iteration  of  the  BOXCAR 

inversion.  The  length  of  the  ticks  indicate  the  relative  importance  of 
each  residual  in  determining  the  solution. 


INLET 


Since  we  have  little  a  priori  information  on  a  structure  model  specific 
to  the  Pahute  Mesa  event  INLET  (11/20/75,  500  kT  estimated),  the  inversion 
procedure  will  follow  the  approach  outlined  in  Test  #  U.  For  the  first  five 
iterations,  we  solve  for  source  and  upper  gradient  parameters.  Using  these 
results  as  a  starting  model,  we  invert  for  source  and  a  two-layered  gradient 
model  with  a  second-order  discontinuity.  The  data  consist  of  vertical 
velocity  waveforms  from  the  three  receivers  listed  in  Table  5»  In  the  first 
five  iterations,  we  hope  to  fit  primarily  travel  time,  so  the  covariance 
matrix  of  the  residuals  is  defined  as  diag(l.O,  1.0,  0.001)  for  waveform, 
amplitude  and  travel  time  residuals  at  S-5  and  S-6,  respectively,  and 
diag(0.1,  0.1,  0.001)  for  the  same  residuals  at  S-7.  The  parameter  weights 
are  defined  as  diag  (0.1,  1.0,  0.01,  1.0)  for  7^,  k,  upper  gradient  and  top 
velocity,  respectively.  In  iterations  6-15,  the  covariance  matrix  of  the 
residuals  is  diag(0.01,  0.1,  0.001)  for  S-5  and  S-6  and  diag(0.1,  0.1,  0.01) 
for  S-7.  The  parameter  weights  are  diag(0.1,  0.1,  0.01,  0.01,  1.0,  0.1), 
where  the  last  two  correspond  to  lower  gradient  and  interface  depth, 
respectively. 

The  results  of  the  inversion  are  listed  in  Table  7,  the  waveforms  are 
plotted  in  Figure  lB,  the  parameter  resolution  matrices  are  shown  in  Figure  19 
and  the  data  Importance  vectors  are  plotted  in  Figure  20.  As  planned,  the 
first  five  iterations  are  controlled  by  the  travel-time  residuals  (Figure  20 ) 
as  well  as  the  amplitude  residual  at  S-7.  The  upper  gradient  parameters  are 
perfectly  resolved  (Figure  19 ),  k  is  fairly  well  resolved,  but  trades  off 
with  k  and  the  upper  gradient.  The  waveforms  (Figure  lR)  change  little  in  the 
first  five  iterations,  except  that  later  arriving  features  are  beginning  to 
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Table  T  -  INLET  Inversion  Results 


Iteration 

Yoo 

( xl010cm3 ) 

k 

(s'1)  (k 

Vo, 

m/s/km) 

(km}s ) 

VOo. 

(km/s/km) 

7j2 

(km) 

RMS  fit 

Start 

3.20 

8.00 

0.83 

2.50 

- 

- 

0.2733 

1 

3.2U 

(0.15) 

U.91 

(1.^9) 

1.21 

(0.42) 

2.64 

(0.08) 

- 

- 

0.4206 

2 

3.3** 

(0.08) 

8.08 

(1.50) 

1.42 

(0.59) 

2.35 

(0.24) 

- 

- 

0.1895 

3 

3.56 
(0.44 ) 

9.60 

(1.65) 

1.41 

(0.42) 

2.4l 

(0.15) 

- 

- 

0.1407 

4 

3.61 

(0.56) 

9.78 

(1.32) 

1.42 

(0.49) 

2.36 

(0.14) 

- 

- 

0.1511 

5 

3.61 

(0.52) 

9-70 

(1.21) 

1.43 

(0.52) 

2.37 

(0.16) 

[0.30] 

[2.50] 

0.1675 

6 

*4.75 

(1.09) 

6.76 

(1.55) 

1.17 

(0.56) 

2.6l 

(0.39) 

1.98 

(0.71) 

3.58 

(0.80) 

0.1781 

7 

5.13 

(0.69) 

7.12 

(0.9M 

1.30 

(0.43) 

2.52 

(0.33) 

2.40 

(0.52) 

3.47 

(0.29) 

0.1521 

8 

5.95 

(1.70) 

6.20 

(1.91) 

1.37 

(0.51) 

2.58 

(0.45) 

3.07 

(0.73) 

2.37 

(2.13) 

0.1269 

9 

6.25 

(0.83) 

6.33 

(0.86) 

1.34 

(0.82) 

2.60 

(0.36) 

2.91 

(0.33) 

3.22 

(2.13) 

0.1491 

10 

6.25 

(0.39) 

6.12 

(1.00) 

1.51 

(0.44) 

2.58 

(0.27) 

3.12 

(0.48) 

3.06 

(0.60) 

0.1369 

11 

6.30 

(0.61) 

6.21 

(1.08) 

1.47 

(0.77) 

2.59 

(0.65) 

3.13 

(0.m6) 

3.06 

(0.03) 

0.1392 

12 

6.57 

(0.56) 

6.64 

(0.95) 

1.27 

(0.74) 

2.58 

(0.42) 

3.31 

(0.16) 

3.18 

(0.52) 

0.1288 

13 

6. U7 
(0.52) 

6.38 

(0.94) 

1.42 

(0.54) 

2.54 

(0.42) 

3.30 

(0.04) 

3.14 

(0.10) 

0.1480 

14 

6.54 

(0.60) 

6.23 

(1.13) 

1.57 

(0.74) 

2.50 

(0.46) 

3.27 

(0.09) 

3.19 

(0.23) 

0.1605 

15 

6.80 

(0.85) 

6.65 

(1.39) 

1.55 

(0.69) 

2.48 

(0.38) 

3.25 

(0.13) 

3.28 

(0.50) 

0.1533 

—  Numbers  in  parentheses  are  the  standard  deviations  of  the  parameter 
changes  that  result  in  the  model  parameters  for  that  iteration. 

—  Numbers  in  square  brackets  are  starting  values  for  iterations  6-15, 
not  results  of  iteration  5. 
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INLET  Inversion  Results 
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Figure  18.  Inversion  results  for  INLET.  Shown  are  synthetic  seismograms  for 
the  starting  model  and  selected  iterations  (above)  compared  with  the 
observed  waveforms  (below).  The  format  is  the  same  as  Figure  15. 
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arise  at  S-7.  The  standard  deviations  of  the  parameter  changes  (Table  7) 


indicate  that  ^  is  not  significantly  changed.  The  decrease  in  top  velocity 
and  increase  in  upper  gradient  msy  indicate  that  near-surface,  high-gradient 
structure  that  is  not  included  in  the  current  parameterization  is  affecting 
the  results.  At  ary  rate,  a  significant  improvement  in  RMS  fit  is  achieved  in 
the  first  five  iterations.  In  iterations  6-15,  upper  gradient  parameters  are, 
once  again,  well  resolved  (Figure  19).  ^  trades  off  with  k,  and  the  lower 
structure  parameters  are  relatively  poorly  resolved.  All  of  the  residuals 
contribute  essentially  equally  to  the  solution  (Figure  20 ),  although  the 
importance  of  the  relatively  simple  waveform  at  S-5  is  somewhat  diminished. 
Features  of  the  synthetic  waveforms  at  S-5  and  S-6  do  not  change  significantly 
in  the  later  iterations,  but  fit  the  observations  reasonably  well  within  the 
inversion  time  window  (Figure  18).  At  S-7,  the  rise  to  the  first  peak  is  well 
modeled,  although  constructive  interference  of  a  second  arrival  broadens  the 
first  peak  of  the  synthetic  beyond  that  observed.  The  interference  of  several 
arrivals  causes  the  trough  to  be  smoothed  into  a  broad,  slowly  rising  feature, 
as  observed  in  the  data. 

The  poor  resolution  of  the  lower  layer  parameters,  the  erratic 
convergence  when  these  parameters  are  added  to  the  inversion,  and  the  final 
values  obtained  for  these  parameters  indicate  either  that  the  three  observed 
waveforms  are  incapable  of  resolving  structure  beneath  about  2.5  km,  or  that 
the  constraint  of  a  second-order  discontinuity  is  too  severe.  At  any  rate,  it 
is  unlikely  that  an  increasing  gradient  below  3  km  is  uniquely  resolved  by 
this  data  set.  The  results  of  iteration  5  indicate  that  the  single  gradient 
model  results  in  source  parameter  estimates  slightly  higher  than  those  of 
Barker,  et  al.  (1985),  which  were  used  as  starting  values.  On  the  other  hand. 
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the  two- layer  nuclei  results  in  an  estimate  of  about  double  that  of  Barker, 
et  al.  (1985),  while  the  k  estimate  is  somewhat  decreased.  Future  inversions 
with  this  data  set  will  free  the  discontinuity  between  layers  to  be  first 
order  in  order  to  determine  whether  the  rise  time  and  amplitude  of  the  first 
peak  are  affected  by  more  reasonable  structures  at  depth. 

MAST 

The  procedure  for  the  inversion  of  waveforms  from  the  Pahute  Mesa  event 
MAST  (6/19/75,  520  kT  estimated)  is  similar  to  that  for  INLET.  The  receivers 
to  be  included  are  listed  in  Table  5,  and  the  covariance  matrix  for  the  first 
five  iterations  is  defined  as  diag(0.01,  0.01,  0.01)  for  the  waveform, 
amplitude  and  travel  time  residuals  of  the  three  stations,  respectively.  The 
parameter  weighting  matrix  is  diag(0.1,  1.0,  0,01,  1.0)  for  k,  upper 
gradient  and  top  velocity,  respectively.  For  iterations  6-15,  the  residual 
covariance  matrix  is  dlag(0.01,  0.1,  0.001)  for  S-5  and  S-6  and  diag(0.01, 
0.01,  0.01)  for  S-7.  The  parameter  weighting  matrix  becomes  diag(0.1,  0.1, 
0.01,  0.01,  0.1,  0.1),  where  the  last  two  correspond  to  lower  gradient  and 
interface  depth,  respectively. 

The  results  of  the  inversion  are  listed  in  Table  B,  and  the  waveforms, 
parameter  resolution  matrices  and  data  importance  vectors  are  plotted  in 
Figures  21-23,  respectively.  As  with  INLET,  the  upper  structure  parameters 
are  well  resolved,  lower  structure  parameters  are  poorly  resolved,  and  u'a,  and. 
k  trade  off  with  k  slightly  better  resolved.  The  interface  depth  is  actually 
better  resolved  than  with  INLET,  but  the  effect  is  to  lower  the  interface  to  a 
point  at  which  the  available  waveforms  are  unlikely  to  resolve  the  structure 


Table  8  -  MAST  Inversion  Results 


Iteration 

too 

(xlO10cm3) 

k  VOL 

(s-1)  (km/s/km) 

(km}s ) 

Vot  p 

(km/s/km) 

Z2 

(km) 

RMS  fit 

Start 

4.70 

7.50 

0.83 

2.50 

- 

- 

0.2646 

1 

4.63 

(0.04) 

6.6  2 
(0.25) 

0.77 

(0.16) 

3.08 

(0.08) 

- 

- 

0.2052 

2 

4.75 

(0.05) 

8.19 

(0.42) 

0.77 

(0.14) 

2.88 

(0.14) 

- 

- 

0.1813 

3 

4.62 

(0.06) 

7.02 

(0.31) 

0.90 

(0.17) 

2.87 

(0.10) 

- 

- 

0.1679 

1* 

4.67 

(0.04 ) 

7.15 

(0.36) 

0.98 

(0.19) 

2.90 

(0.11) 

- 

- 

0.1615 

5 

4.66 

(0.06) 

7.02 

(0.37) 

1.12 

(0.25) 

2.91 

(0.13) 

[0.20] 

[2.00 1 

0.1603 

6 

4.69 

(0.20) 

5.89 

(1.19) 

1.02 

(0.44) 

3.07 

(0.58) 

0.35 

(0.34) 

2.64 

(0.41) 

0.2197 

7 

5.01 

(0.57) 

6.58 

(1.48) 

0.99 

(0.79) 

2.96 

(0.85) 

0.35 

(0.32) 

2.75 

(0.6l) 

O.I565 

8 

5-03 

(0.31) 

6.60 

(0.52) 

0.99 

(0.19) 

2.95 

(0.18) 

0.35 

(0.01) 

2.77 

(0.59) 

0.1347 

9 

4.95 

(0.15) 

6.34 

(0.45) 

1.09 

(0.20) 

2.95 

(0.18) 

0.39 

(0.05) 

2.79 

(0.18) 

O.I691 

10 

5.10 

(0.24) 

6.36 

(0.46) 

1.12 

(0.33) 

2.93 

(0.16) 

0.42 

(0.05) 

3.21 

(0.88) 

0.1527 

11 

5.18 

(0.65) 

5.81 

(1.46) 

1.18 

(0.30) 

2.91 

(0.29) 

0.55 

(0.51) 

3.63 

(1.76) 

0.2067 

12 

5.33 

(0.27) 

5.94 

(0.44) 

1.26 

(0.27) 

2.80 

(0.38) 

O.56 

(0.01) 

3.62 

(0.01) 

0.1624 

13 

5.78 

(0.60) 

5.81 

(0.67) 

1.36 

(0.44) 

2.80 

(0.28) 

0.63 

(0.11) 

4.44 

(1.31) 

0.1517 

lU 

6.00 

(0.42) 

6.02 

(0.47) 

1.25 

(0.53) 

2.89 

(0.52) 

0.63 

(0.01) 

4.35 

(0.23) 

0.1213 

15 

6.41 

(0.96) 

5.80 

(0.43) 

1.48 

(0.84) 

2.70 

(0.59) 

0.60 

(0.06) 

4.51 

(0.26) 

0.1074 

—  Numbers  in  parentheses  are  the  standard  deviations  of  the  parameter 
changes  that  result  in  the  model  parameters  for  that  iteration. 

—  Numbers  in  square  brackets  are  starting  values  for  iterations  6-15, 
not  results  of  iteration  5. 
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Figure  21.  Inversion  results  for  MAST.  Shown  are  synthetic  seismograms  for 
the  starting  model  and  selected  iterations  (above)  compared  with  the 
observed  waveforms  (below).  The  format  is  the  same  as  Figure  15. 


The  parameter  resolution  matrices  for  each  iteration  of  the  inversion  of  MAST  waveforms 
gradient  parameters  were  not  considered  in  the  first  five  iterations. 


beneath  it  (see  Table  8).  Also  unlike  INLET,  the  lower  gradient  remains 
smaller  than  the  upper  gradient,  suggesting  no  need  for  a  first-order 
discontinuity  at  depth.  Except  for  the  improvement  in  RMS  fit  in  the  last 
iteration,  a  single  gradient  model  would  appear  as  reasonable  as  the  two-layer 
model.  The  last  iteration  yields  improved  fits  to  the  amplitudes  of  the  first 
peaks  (Figure  21),  begins  to  account  for  the  double  arrival  in  the  first  peak 
at  S-T,  and  includes  multiple  interference  features  in  the  first  trough  at  S-6 
and  S-T.  The  data  importance  vectors  (Figure  23)  suggest  that  the  simple 
waveform  at  S-5  is  relatively  unimportant  to  the  solution.  The  Importance  of 
the  other  residuals  follows  the  weighting  scheme. 

If  we  consider  that  the  single  gradient  model  is  sufficient,  the  results 
of  iteration  5  (Table  8)  indicate  that  there  is  no  significant  change  in 
from  the  value  obtained  by  Barker,  et  al.  (1985)  (the  starting  model).  The 
value  obtained  for  k  is  slightly  decreased.  Once  again,  the  upper  gradient  is 
increased  and  the  top  velocity  is  decreased  relative  to  the  Hartzell,  et  al. 
(1983)  model,  although  this  could  be  an  effect  of  ignoring  near-surface, 
high-gradient  structure.  With  the  inclusion  of  the  lower  layer  (whether 
actually  resolved  or  not),  the  estimate  of  is  increased,  that  of  k  is 
decreased,  and  the  value  of  the  upper  gradient  is  further  increased.  The  top 
velocity  is  not  changed  significantly  from  the  starting  model. 
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SCOTCH 


For  the  inversion  of  the  Pahute  Mesa  event  SCOTCH  (5/23/67,  155  kT 
announced),  there  are  only  two  local  observations  outside  the  spall  zone  (see 
Table  5)«  Once  again,  following  the  procedure  outlined  above,  we  invert  first 
for  source  and  upper  structure  parameters,  setting  the  parameter  weighting 
matrix  to  be  diag(0.1,  1.0,  0.01,  1.0)  for  ^co,  k,  upper  gradient  and  top 
velocity,  respectively.  The  covariance  matrix  for  both  receivers  is 
diag(0.01,  0.01,  0.01)  for  waveform,  amplitude  and  travel-time  residuals, 
respectively.  For  iterations  6-15,  we  include  the  parameters  of  the  lower 
layer,  and  define  the  weighting  matrix  as  diag(0.1,  1.0,  0.01,  0.01,  1.0, 

0.1 ),  where  the  last  two  correspond  to  lower  gradient  and  interface  depth, 
respectively.  The  residual  covariance  matrix  is  adjusted  to  diag(0.01,  0.1, 
0.001 )  for  both  receivers. 

The  inversion  results  are  listed  in  Table  9  and  the  waveforms,  parameter 
resolution  matrices  and  data  importance  vectors  are  plotted  in  Figures  2U-26. 
The  parameter  resolution  matrix  (Figure  25)  indicates  that,  as  before,  the  top 
velocity  is  perfectly  resolved  throughout  the  inversion.  The  upper  gradient 
is  well  resolved  except  in  iterations  3-5,  in  which  it  is  hardly  resolved  at 
all.  As  usual,  and  k  trade  off,  this  time  about  equally.  The  lower 
structure  parameters  are  somewhat  poorly  resolved,  although  in  a  few 
iterations,  the  interface  depth  is  resolved.  Since  there  are  only  two 
receivers,  each  of  the  residuals  is  equally  important  in  defining  the  solution 
(Figure  26). 

Within  the  first  two  or  three  iterations,  the  upper  gradient  and  top 
velocity  have  essentially  reached  their  final  values  (see  Table  9).  These 
parameters  are  well  resolved  by  the  data.  Hy  the  fifth  iteration,  a  strong  pP 


Table  9  -  SCOTCH  Inversion  Results 


Iteration 

H'co 

( xlO10cra^) 

k 

(s'1)  ( 

[km/s}km)  (km}s) 

Vap 

(km/s/km) 

Z2 

1  (km) 

RMS  fit 

Start 

1.30 

12.00 

0.83 

2.50 

- 

- 

0.2624 

1 

1.45 

(0.22) 

10.85 

(0.81*) 

1.22 

(0.39) 

2.35 

(0.14) 

- 

- 

0.2412 

2 

2.22 

(0.25) 

12.1*0 

(0.28) 

1.76 

(0.1*5) 

2.01 

(0.17) 

- 

- 

0.1849 

3 

2.22 

(0.19) 

12.1*0 

(0.26) 

1.75 

(0.02) 

1.83 

(0.08) 

- 

- 

0.1592 

U 

2.23 

(0.18) 

12.1*1+ 

(0.28) 

1.75 

'0.05) 

1.84 

(0.06) 

- 

- 

0.1601 

5 

2.2** 

(0.15) 

12.1*5 

(0.25) 

1.75 

(0.05) 

1.8l 

(0.06) 

[0.20] 

12.00] 

0.1794 

6 

1.76 

(0.68) 

11.1*6 

(0.86) 

1.56 

(0.23) 

2.01 

(0.18) 

2.21 

(0.64) 

2.48 

(0.49) 

0.2707* 

T 

1.67 

(0.U2) 

11.79 

(0.71) 

1.89 

(0.35) 

1.92 

(0.37) 

3.94 

(0.66) 

2.03 

(0.41) 

O.2365 

8 

1.63 

(0.1*5) 

11.1*6 

(1.19) 

1.85 

(0.52) 

1.94 

(0.28) 

3.91 

(0.16) 

2.37 

(0.27) 

0.2581 

9 

2.1*7 

(0.81*) 

13.60 

(0.97) 

1.97 

(0.44) 

1.8l 

(0.25) 

3.85 

(0.14) 

2.53 

(0.18) 

0.3326 

in 

2.68 

(0.69) 

II.85 

(0.75) 

1.84 

(0.29) 

1.84 

(0.34) 

6.74 

(1.97) 

3.22 

(0.64) 

0.1692 

ll 

2.52 

(0.62) 

12.63 

(0.68) 

1.83 

(0.47) 

1.88 

(0.35) 

6.68 

(0.18) 

3.65 

(0.99) 

0.1953 

12 

2.28 

(0.43) 

12.65 

(0.95) 

2.07 

(0.34) 

1.78 

(0.21) 

7.07 

(0.25) 

3.67 

(0.04) 

0.1819 

13 

2.17 

(0.54) 

10.88 

(1.1*6) 

1.86 

(0.33) 

1.96 

(0.36) 

7.71 

(0.34) 

3.64 

(0.05) 

0.2159 

lU 

1.1*3 

(1.06) 

13.10 

(0.99) 

1.88 

(0.37) 

1.90 

(0.24) 

8.50 

(0.31) 

2.30 

(0.90) 

0.2759 

15 

2.78 

(0.86) 

8.90 

(1.24) 

1.84 

(0.20) 

1.90 

(0.06) 

9.23 

(0.19) 

3.53 

(O.58) 

0.2049 

—  Numbers  in  parentheses  are  the  standard  deviations  of  the  parameter 
changes  that  result  in  the  model  parameters  for  that  iteration. 

—  Numbers  in  square  brackets  are  starting  values  for  iterations  6-15, 
not  results  of  iteration  5. 

The  time  window  of  S-4  was  lengthened  at  iteration  6. 
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SCOTCH  Inversion  Results 


4.13  1cm 
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Figure  2U.  Inversion  results  for  SCOTCH.  Shown  are  synthetic  seismograms  for 
the  starting  model  and  selected  iterations  (above)  compared  with  the 
observed  waveforms  (below).  The  format  is  the  same  as  Figure  15* 
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phase  may  be  identified  in  the  synthetics  for  station  S-3A  (Figure  2*0  which 
correlates  well  with  an  arrival  in  the  observed  waveform  Just  outside  of  the 
inversions  time  window.  Similarly,  the  high  upper  gradient  allows  the  upgoing 
direct  wave  at  station  S-4  to  separate  in  time  from  the  first-arriving  diving 
ray,  and  is  apparent  as  the  second  peak  in  the  observed  waveform.  Because 
this  arrival  is  well  modeled,  the  inversion  time  window  for  station  S-U,  which 
originally  included  only  the  first  peak  and  trough,  was  lengthened  to  that 
indicated  in  Figure  2U.  (This  also  accounts  for  the  increase  in  RMS  fit  at 
iteration  6  in  Table  9).  The  pP  ray  at  S-U  is  quite  strong  in  the  synthetics 
Just  outside  the  inversion  time  window,  but  is  less  well  correlated  with  the 
observed  waveform.  From  the  sixth  iteration,  the  depth  to  the  lower  layer 
increases  to  about  3.5  km,  while  the  lower  gradient  continuously  increases. 
This  is,  once  again,  indicative  that  the  constraint  of  a  second-order 
discontinuity  may  be  in  error,  and  that  a  first-order  velocity  discontinuity 
beneath  the  source  may  be  required.  Throughout  the  inversion,  and  k  vary 
notably,  reducing  confidence  in  the  value  obtained  at  ary  particular 
iteration.  Although  the  lower  gradient  has  not  yet  converged  to  a  final  value 
(and  may  not)  by  the  fifteenth  iteration,  the  model  obtained  after  this 
iteration  provides  a  good  fit  to  the  first  peak  amplitudes  (Figure  2h)  and  a 
reasonably  improved  RMS  fit  (Table  9)*  The  value  of  obtained  for  this 
iteration  is  more  than  double  the  value  obtained  by  Barker,  et  al.  (1985)  (the 
starting  model).  The  corresponding  value  of  k  is  the  smallest  of  any  of  the 
iterations,  and  is  significantly  smaller  than  the  value  of  Barker,  et  al. 
(1985). 
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SUMMARY  AND  CONCLUSIONS 


We  have  introduced  a  method  for  the  simultaneous  inversion  of  near-field 
waveforms  for  source  and  structure  parameters.  For  this  report,  the  source 
was  represented  by  the  Helmberger-Hadley  (1981)  RDP  with  B  held  constant. 
However,  any  of  a  number  of  other  effective  source  functions  could  be 
substituted.  The  structure  model  is  parameterized  as  a  layered  stack,  in 
which  the  free  parameters  are  the  velocity  at  the  top  of  each  layer,  the 
velocity  gradient  within  the  layer  and  the  depth  to  the  top  of  the  layer.  For 
each  observation,  the  residuals  consist  of  the  normalized  cross-correlation 
coefficient  between  the  data  and  synthetics  (a  measure  of  waveform  fit),  the 
difference  of  the  normalizations  (a  measure  of  absolute  amplitude),  and  the 
time  lag  to  the  peak  of  the  cross-correlation  (a  measure  of  absolute  travel 
time).  Numerical  partial  derivatives  are  computed  using  asymptotic 
generalized  ray  synthetics,  and  the  problem  is  solved  using  an  iterative 
generalized  inverse. 

Test  inversions  have  been  performed  for  simple  structure  models.  The 
inversion  for  half space  velocity  is  quite  robust,  but  an  inversion  for  a 
single  gradient  indicates  that  the  nonlinearity  of  the  problem  can  cause  slow 
convergence,  or  even  temporary  divergence,  of  the  solution.  Methods  of 
nonlinear  optimization  could  minimize  this  problem,  but  have  not,  as  yet,  been 
implemented.  In  solving  for  both  surface  velocity  and  gradient,  the 
parameters  trade  off  and  the  inversion  is  unstable.  However,  truncating 
singular  values  to  damp  poorly  resolved  parameter  changes,  and  preferentially 
weighting  the  parameters  and  residuals,  allows  the  inversion  to  slowly 
approach  the  correct  solution.  In  a  final  test  inversion  for  the  source  and  a 


two-gradient  structure  model,  the  inversion  is  far  more  sensitive  to  the 
parameters  of  the  upper  gradient  layer  than  to  either  the  source  or  lower 
layer  parameters.  A  procedure  was  outlined  in  which  we  solve  in  early 
iterations  for  estimates  of  source  and  upper  gradient  parameters  only,  then 
use  these  as  starting  models  in  further  iterations  for  the  total  model.  Even 
with  this  approach,  however,  the  lower  gradient  parameters  could  be  resolved 
only  when  the  range  of  "observed"  data  was  extended. 

The  inversion  was  applied  to  near-field  waveforms  recorded  for  four 
Pahute  Mesa  explosions  to  obtain,  for  each,  estimates  of  the  source  and  of  a 
two-gradient  structure  model  with  a  second-order  discontinuity.  The  resulting 
structure  models  are  summarized  in  Figure  27*  For  BOXCAR,  the  structure  model 
was  not  significantly  different  from  the  starting  model,  which  was  based  on 
the  Pahute  Mesa  model  of  Hartzell,  et  al.  (1983).  The  structure  models 
obtained  for  INLET,  MAST  and  SCOTCH,  however,  include  a  higher  gradient  in  the 
upper  layer,  and  poor  resolution  of  the  lower  layer.  In  fact,  the  lower 
gradient  for  INLET  and  SCOTCH  is  actually  greater  than  the  upper  gradient, 
suggesting  that  the  constraint  that  the  interface  have  only  a  second-order 
discontinuity  may  be  incorrect.  Future  inversions  of  these  data  sets  will 
free  the  model  to  have  a  first-order  velocity  discontinuity  at  depth. 

Although  these  results  should  be  considered  preliminary,  the  source  parameters 
obtained  for  the  four  events  all  give  estimates  of  4’^  higher  than  those 
determined  by  Barker,  et  al.  (1985)  using  the  Hartzell,  et  al.  (1983) 
structure  model.  As  expected  for  increased  estimates  of  4*  the  estimates  of 
k  are  decreased  somewhat  from  the  Barker,  et  al.  (1985)  results  for  INLET, 

MAST  and  SCOTCH.  However,  the  estimate  of  k  for  BOXCAR  is  increased. 

Future  improvements  in  the  simultaneous  inversion  procedure  might  include 
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Figure  27.  A  comparison  of  velocity  structure  models  resulting  from  the 

waveform  inversions  of  the  Pahute  Mesa  events  BOXCAR,  INLET,  MART  and 
SCOTCH.  The  model  for  BOXCAR  is  very  similar  to  that  of  Hartzell,  et  al. 
(19B3).  The  other  models  suggest  higher  velocity  gradients  in  the  upper 
3-1*  km.  The  deeper  portions  of  these  models  are  poorly  resolved  and  may 
reflect  the  constraints  placed  on  the  parameters. 
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the  implementation  of  Shaw  and  Orcutt’s  (1985)  "jumping"  method,  as  well  as 
nonlinear  optimization  techniques  such  as  the  "conjugate  gradient"  method. 
Changes  in  parameterization  include  the  specification  of  surface  velocity  and 
water  table  depth,  when  known,  for  Pahute  Mesa  events,  as  well  as  allowing 
first-order  velocity  discontinuities  at  depth.  Using  a  quadratic 
parameterization  of  the  explosion  RDP  (e.g.  Mueller  and  Murphy,  1971)  will 
require  inclusion  of  depth-dependent  attenuation  in  the  structure  model.  It 
will  be  of  interest,  however,  to  determine  how  the  choice  of  source  function 
affects  the  structure  model  obtained,  and  to  investigate  the  relative  waveform 
matches  of  the  two  approaches.  Applications  to  data  from  other  sites  might 
include  surface  and  borehole  records  of  Yucca  Flats  events  modeled  by  Wallace 
and  Barker  (1986),  as  well  as  near-field  data  from  non-NTS  shots  such  as 
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Part  II 

Q  AND  YIELD  ESTIMATES  FROM  THE 
PAHUTE  MESA  TEST  SITE 


R.  W.  Burger,  T.  Lay  and  L.  J.  Burdick 


ABSTRACT 


Attenuation  models,  with  and  without  frequency  dependence,  have  been 
developed  through  analysis  of  time-domain  amplitude  measurements  and 
teleseismic  spectral  shape  data  from  Pahute  Mesa  nuclear  explosions.  The 
time-domain  analysis  is  based  on  a  near-field  to  far-field  amplitude 
comparison.  The  near-field  amplitude  information  is  incorporated  in  two 
parameterized  explosion  source  models  (Mueller-Murphy  and  Helmberger-Hadley ) 
based  on  analyses  of  near-field  data.  The  teleseismic  amplitude  observations 
are  from  a  large  dataset  of  WWSSN  short-period  analog  recordings.  For  the 
narrow-band  time-domain  data,  the  various  source  and  attenuation  models  are 
indistinguishable.  We  utilize  the  spectral  shape  data  in  the  0.5  to  U  Hz  band 
as  a  constraint  on  the  source-attenuation  models  at  higher  frequencies, 
concluding  that  either  source  model,  when  convolved  with  the  appropriate 
frequency  dependent  Q  model,  can  be  consistent  with  both  the  near-field  and 
far-field  time-domain  amplitudes  and  the  spectral  shape  data.  Given  the 
trade-off  between  source  and  attenuation  models  and  the  similarity  of  the 
different  source  models  in  the  0.5  to  U  Hz  band,  it  is  difficult  to  clearly 
prefer  one  source  model  over  the  other.  The  Mueller-Murphy  model  is  more 
consistent  with  surface  wave  amplitude  measurements,  because  of  larger 

predicted  long-period  energy  levels.  Whether  or  not  frequency  dependence  is 

* 

included  in  the  attenuation  model,  the  value  of  t  near  1  Hz  is  about  1.0  sec 
assuming  the  Mueller-Murphy  source,  or  0.8  sec  assuming  the  Helmberger-Hadley 
source  model.  This  0.2  sec  difference  results  from  greater  1  Hz  energy  levels 
for  the  Mueller-Murphy  source  model.  Adopting  an  average  attenuation  model, 
predicted  amplitudes  and  yields  are  shown  to  be  within  the  uncertainty  of  the 
data  for  all  the  events  analyzed. 


INTRODUCTION 


Understanding  the  nature  and  effect  of  the  earth's  anelasticity  on 
seismic  booty-  waves  is  an  important  problem  in  seismic  yield  estimation. 
Reliable  estimates  of  yield  based  on  network -averaged  short-period  m^  require 
detailed  knowledge  of  P-wave  attenuation  in  the  earth.  Furthermore,  it  is  now 
recognized  that  the  earth's  anelasticity  has  significant  lateral  variations. 

ft 

Thus  it  is  necessary  to  estimate  the  deviation  of  t  from  the  global  average 
in  the  short-period  body  wave  band  at  a  given  test  site  to  predict  the 
m^ -yield  relation  baseline  for  that  site.  Bache  et  al.  (I9B5,  1986)  and  Der 
et  al.  (1985)  have  used  spectral  shapes  and  decay  rates  of  teleseismic  P  waves 
from  explosions  to  estimate  the  level  and  frequency  dependence  of  Q.  However, 
spectral  methods  for  estimating  Q  have  several  potential  sources  of  error 
(Cormier,  1982).  Most  important  is  that  spectral  decay  methods  require  a 
priori  knowledge  of  the  source  spectrum.  Errors  in  the  assumed  source 
spectral  falloff  rate  translate  into  a  bias  in  the  attenuation  measurement. 
This  presents  a  problem  in  the  spectral  decay  studies  of  explosion  signals 
given  that  explosion  source  models  with  different  spectral  falloff 
characteristics  have  been  proposed  at  various  times  (e.g.  Haskell,  19&7; 
Mueller  and  Murphy,  1971;  von  Seggern  and  Blandford,  1972;  and  Helmberger 
and  Hadley,  1981) 

A  common  misconception  is  that  these  source  models  can  be  effectively 
characterized  by  their  high  frequency  spectral  falloff  rates  (f-2,f-^,  or 
f-1*).  However,  we  show  in  this  work  that  this  is  not  necessarily  the  case. 

For  the  wide  range  of  yields  investigated  in  this  study  (155  to  1300  Kt*,  the 
Helmberger  and  Hadley  (1981)  source  which  has  an  asymptotic  spectral  decay 
rate  of  f~  ,  has  spectral  properties  in  the  band  of  interest  (0.5-**  Hz) 


similar  to  the  Mueller  and  Murphy  (1971)  source,  which  has  an  f“^  asymptotic 
decay.  We  also  show  that  these  source  nndels,  when  convolved  with  an 
appropriate  attenuation  model,  are  indistinguishable  in  WWSSN  short-period 
data.  It  is  likely  that  none  of  these  source  models  is  completely  accurate. 
The  spectrum  may  initially  decay  at  one  rate  and  then,  at  a  higher  frequency, 
begin  to  decay  at  a  faster  rate. 

Time-domain  modeling  of  explosion  signals  provides  an  alternate  method 
to  the  spectral  shape  procedures  for  estimating  attenuation.  However, 
separate  analysis  of  time-domain  and  spectral  decay  information  may  produce 
incompatible  results,  given  the  different  intrinsic  sensitivity  of  each 
procedure  (Cormier,  I9S2).  In  this  investigation,  we  have  measured  the 
attenuation  of  P  waves  from  explosions  at  the  Pahute  Mesa  test  site  utilizing 
time-domain  amplitude  data  in  the  0.5-2  Hz  band  (WWSSN  short-period)  and 
spectral  decay  information  as  a  constraint  on  the  source-attenuation  models  at 
frequencies  up  to  U  Hz. 

To  eliminate  the  trade-off  between  source  and  attenuation,  several 
investigators  have  used  source  canceling  methods  to  isolate  the  effects  of 
propagation  through  the  anelastic  earth.  Jordan  and  Sipkin  (1977)  and  Sipkin 
and  Jordan  (1979 ,  1980)  among  others  have  used  multiple  ScS  waves  as  a  direct 
measure  of  the  average  shear  wave  Q  in  the  mantle.  Since  each  leg  of  the  ScS 
arrival  should  have  the  same  source  excitation,  differences  in  the  frequency 
content  of  successive  multiples  can  be  attributed  to  attenuation.  Burdick 
(1985)  used  the  relative  frequency  content  of  ScS  and  ScP  waves  to  estimate 
the  average  mantle  Q.  These,  and  other  phase  pairs  such  as  sP-sS  (Burdick, 
1978)  and  P-PP  (Shore,  1983),  provide  an  unbiased  estimate  of  attenuation.  A 
recent  review  on  the  subject  of  bocfy-wave  attenuation  was  presented  by  Cormier 
(1982). 
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Phase  pairs  for  which  source  canceling  studies  like  those  described 
above  might  be  carried  out  are  not  commonly  observed  from  explosions. 

However,  nuclear  sources  have  been  well  recorded  in  the  near-field  and 
near-regional  field.  This  allows  an  alternate  procedure  of  determining  the 
initial  source  excitation  from  the  near-field  data  and  estimating  Q  by 
matching  the  observed  teleseismic  amplitudes.  This  type  of  investigation  has 
been  made  previously  by  Frasier  and  Filson  (1972),  Helmberger  and  Hadley 
(1981),  and  Burdick  et  al.  (1984). 

Near-field  records  from  explosions  at  Pahute  Mesa  have  been  studied 
extensively.  Yield-scaled  source  models  have  been  developed  from  this  data 
using  two  different  approaches.  Mueller  and  Murphy  (1971)  compared  yield 
scaling  exponents  as  a  function  of  frequency  for  events  of  various  sizes 
recorded  at  a  common  station.  Murphy  (1977)  later  refined  their 
representation  into  what  is  referred  to  in  this  paper  as  the  Mueller-Murply 
source  model.  Helmberger  and  Hadley  (1981)  developed  their  f~^  model  by 
forward  modeling  of  a  small  set  of  observations,  but  they  did  not  present 
yield  scaling  relations.  Pahute  Mesa  yield  scaling  relations  for  the  latter 
source  representation  were  later  developed  by  Barker  et  al.  (1985)  along  with 
source  models  for  a  wider  range  of  events.  These  two  different  types  of 
studies  account  for  the  complex  process  of  near-field  wave  propagation  in  very 
different  ways,  and  assume  different  parametric  forms  for  the  source  time 
history  (spectrum).  We  consider  both  the  Mueller-Murphy  (1971)  and 
Helmberger-Hadley  (1981)  source  models  in  our  modeling  to  illustrate  the 
trade-offs  between  source  and  attenuation  models  in  the  0.5  to  4  Hz  band,  and 
to  appraise  the  various  claims  in  the  literature  that  one  model  is  superior  to 
the  other. 

It  is  widely  recognized  that  the  effective  value  of  Q  for  the  mantle 


varies  with  frequency  within  the  body  wave  band  (0.2  to  20  sec).  However,  it 
is  common  practice  to  assume  that  within  the  frequency  band  of  any  one  type  of 
narrow-band  data,  Q  is  independent  of  frequency.  The  variation  of  Q  with 
frequency  is  inferred  by  comparing  the  results  of  several  different  types  of 
studies  (e.g.  Burdick,  1985).  In  this  work,  we  have  considered  both  frequency 
independent  (Futterraan,  19^2)  and  frequency  dependent  models  for  attenuation 
filters.  The  frequency  dependent  attenuation  model  that  we  use  is  the 
standard  linear  solid  model  described  by  Minster  (1978a,  1978b).  Our  goal  was 
to  find  whether  either  source  model  or  either  attenuation  model  did  a  superior 
Job  of  fitting  the  combined  teleseismic  spectral  and  amplitude  data  sets  from 
Pahute  Mesa.  Our  conclusion  is  that  any  of  the  four  source-attenuation  model 
pairs  can  satisfy  short-period  WWSSN  observations.  But,  when  using  spectral 
shape  data  in  the  0.5  to  ^  Hz  band  as  an  additional  constraint,  we  find  it 
necessary  to  include  a  frequency  dependent  attenuation  model.  Furthermore; 
for  the  Minster  Q  models,  both  source  models  are  consistent  with  the  amplitude 
and  spectral  shape  data  in  this  band. 


MODELING  APPROACH 

The  basic  approach  taken  in  this  investigation  is  to  attempt  to 
simultaneously  match  observed  teleseismic  amplitudes  and  magnitudes  by 
computing  synthetic  seismograms  using  source  models  which  fit  near-field  data. 
The  parameters  that  we  adjust  to  fit  the  data  are  the  free  parameters  in  two 
standard  attenuation  models.  The  reason  for  matching  both  the 
period-corrected  amplitude  (magnitude)  and  the  uncorrected  amplitude  is  to 
obtain  a  reasonable  measure  of  how  well  the  source  and  attenuation  models 
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predict  the  average  period  of  the  observed  P  waves.  We  also  ensure  that  the 
synthetic  waveforms  agree  with  the  data. 

Observations.  The  observations  that  we  will  attempt  to  match  are  of  two 
different  types.  The  first  type  are  teleseismic  short-period  P-wave  amplitude 
and  magnitude  observations  as  measured  from  WWSSN  film  chips.  A  large  data 
base  of  such  measurements  from  25  Pahute  events  was  processed  in  a  thorough 
and  uniform  fashion  by  Lay  (1985).  He  measured  over  1200  amplitude  values  and 
obtained  from  them  average  log  amplitudes  and  average  magnitudes.  The 
amplitude  and  magnitude  values  to  be  modeled  here  are  given  in  Table  1.  The 
amplitude  that  we  model  is  that  of  the  first  upswing  to  first  downswing  of  the 
explosion  waveform.  This  amplitude  is  not  influenced  by  the  pP  arrival  for 
these  events  (Lay,  1985 ),  but  the  pP  arrival  may  influence  the  period 
measurement.  Thus  the  synthetic  magnitudes  may  be  somewhat  affected  by 
incorrect  pP  delay  time  estimates.  The  band  of  WWSSN  short-period  data  is  0.5 
to  2  Hz.  Therefore  the  time-domain  amplitude  data  discussed  in  this  paper 
(and  apy  possible  spectral  data  from  these  instruments)  are  sensitive  only  to 
the  spectral  level  in  this  band. 

The  second  type  of  observations  that  we  wish  to  match  are  spectral  shape 
observations  in  the  0.5-**  Hz  band.  Unfortunately,  body  wave  spectra  which  are 
hand  digitized  from  analog  records  are  not  reliable  at  high  frequencies,  so  we 
cannot  use  the  spectra  of  the  P  waves  studied  by  Lay  (1985).  Der  et  al. 

(1985)  analyzed  the  spectral  shapes  of  Pahute  Mesa  events  using  several 
different  broad-band  recording  systems.  They  assumed  the  von  Seggern  and 

Blandford  (1972)  source  function  and  determined  a  frequency  independent  value 

£ 

of  t  in  the  0.5-**  Hz  band  of  0.1*  to  0.5  sec  for  Pahute  Mesa.  We  have 

computed  a  synthetic  spectrum  using  this  information  and  we  here  treat  it  as 
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Table  1 


PAHUTE  MESA  EVENTS  STUDIED  * 


Yield 

K 

Y 

Event 

Date 

Log  Amp 

"‘b 

(Kt) 

(s-1) 

(xlO10  , 

ALMENDRO 

6/66/73 

2.U8+0.09 

5.96+0.10 

7.5 

5.7 

BENHAM 

12/19/68 

2.62+0.09 

6. 16+0.09 

1150 

BOXCAR 

L/26/68 

2.59+0.09 

6.12+0.11 

1300 

6.5 

12.0 

GREELEY 

12/20/66 

2.5^+O.K) 

6.09+0.11 

870 

HALFBEAK 

6 /30/66 

2.31+0.09 

5.81+0.09 

365 

9.0 

3.8 

INLET 

11/20/75 

2.23+0.11 

5.71+0.11 

8.0 

3.2 

MAST 

6/19/75 

2.31+0.11 

5.81+0.12 

7-5 

4.7 

SCOTCH 

5/23/67 

1.93+0.10 

5.39+0.11 

155 

12.0 

1.3 

*  The  ampl 

itudes  and 

magnitudes  are 

from  Lay 

(1985). 

The 

Helmberger-Hadley  source  parameters  K  and  are  from 
Barker  et  al.  (1985). 
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an  observation.  We  will  also  consider  the  results  from  the  spectral  studies 
of  Bache  et  al.  (1986).  Th^y  obtained  Pahute  Mesa  spectra  and  magnitude  data 
from  the  UKAEA  and  NORSAR  arrays  and  modeled  it  using  the  Mueller  and  Murphy 

(1971)  source  model.  They  solved  for  frequency  independent  and  frequency 

* 

dependent  values  of  tQ,  concluding  that  the  data  requires  a  frequency 
dependent  Q  operator.  The  spectral  shape  data  from  these  two  studies  are 
quite  consistent  with  each  other,  but  there  is  a  general  difference  in  the 
interpretation  (Bache  et  al. ,  1986).  We  will  constrain  the  models  developed 
here  from  the  time-domain  amplitude  data  to  match  the  spectral  shape  data  as 
well. 


Attenuation  Models.  Two  of  the  most  commonly  used  attenuation  models  in 

bocfy  wave  seismology  are  the  frequency  independent  Futterman  (1962)  and  the 

frequency  dependent  standard  linear  solid  (Minster,  1978a,  b)  models.  The 

Futterman  model  is  parameterized  by  a  single  value,  t*,  obtained  by 

integrating  the  slowness  divided  by  Q  over  the  raypath.  The  standard  linear 

solid  attenuation  model  is  parameterized  by  three  values.  These  are  the 

minimum  value  of  Q,  the  high  frequency  characteristic  relaxation  time,x  ,  and 

m 

the  low  frequency  relaxation  tirae,T  •  T  controls  the  increase  of  Q  at  short 

M  m 

periods  and  x  controls  the  increase  at  long  periods.  Since  we  will  be 
M 

modeling  short-period  observations,  we  simply  fix  TM  at  1000  sec  for  all 

calculations.  The  attenuation  model  is  then  reparameterized  in  terms  of  a 

•  * 

long  period  value  of  t  ,  designated  by  tQ,  and  T^. 

A  reasonable  range  for  t*  can  be  established  from  normal  mode  and 

£ 

long-period  multiple  ScS  data.  Anderson  and  Given  (1982)  give  values  of  tQ 

based  on  normal  mode  data  that  range  from  0.9  to  1.2  sec  between  distances  of 

30°  and  SO0  at  a  period  of  100  sec.  Sipkin  and  Jordan  (1979)  found  that  good 
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quality  multiple  ScS  data  gives  slightly  lower  values  of  Q  than  the  mode  data. 
Burdick  (1982)  found  that  their  observed  average  value  for  whole  mantle  Q 

ft 

translates  to  a  t  of  about  1.3  sec.  Thus  a  reasonable  range  of  values  to 
consider  is  0.9  to  1.3  sec. 

Source  Models.  The  explosion  source  models  which  we  will  use  in  our 
study  are  the  Mueller -Mur pity  tuff  source  and  the  Helmberger-Hadley  source 
representations.  MurpVy  (1977)  provides  the  yield  scaling  relations 
appropriate  for  Pahute  Mesa.  The  form  of  these  relations  is  based  on 
theoretical  and  observational  considerations.  Th^y  provide  an  analytic 
expression  for  the  source  as  a  function  of  yield,  burial  depth  and  material 
properties.  We  assume  that  the  material  properties  are  relatively  constant  at 


Pahute  Mesa,  so  that  as  long  as  yield  is  known  we  can  derive  the  source  time 
history  and  its  spectrum.  Five  events  in  the  Pahute  Mesa  data  set  had 
announced  yields-,  SCOTCH,  HALFBEAK,  GREELEY,  BOXCAR  and  BENrIAM.  The 
announced  yields  are  given  in  Table  1.  In  the  following,  we  use  the  seismic 
observations  from  these  five  events  to  obtain  average  Futterman  and  Minster  Q 
operators.  We  then  predict  yields  for  the  remainder  of  the  events  assuming 
that  these  models  are  correct. 

The  Helmberger-Hadley  source  model  is  based  on  forward  modeling  studies 
of  near-field  data  from  Pahute  Mesa  events  (Hartzell  et  al. ,  1983;  Barker  et 
al.  ,  1985).  The  band  of  their  data  is  1  to  5  Hz,  so  that  their  source  models 
are  appropriate  to  use  in  our  investigation.  Yield  scaling  laws  were 
developed,  but  are  empirical  in  nature  and  based  on  few  data  points.  It  is 
thus  roost  reasonable  to  use  the  source  models  developed  for  particular  events 
whenever  possible  and  we  shall  do  so  in  the  following.  Six  events  were 
forward  modeled  with  good  success  in  these  studies.  They  were  BOXCAR,  SCOTCH, 
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ALMENDRO,  HALFBEAK,  INLET  and  MAST.  The  Helmberger-Hadley  source  parameters  K 
(corner  frequency)  and  ¥  (DC  level)  are  given  in  Table  1.  Again,  we  use 
these  six  events  to  establish  the  attenuation  operators  and  then  predict  the 
yields  of  the  remaining  events. 

Three  events  with  announced  yields  have  been  forward  modeled ;  BOXCAR, 
SCOTCH  and  HALFBEAK.  We  concentrate  primarily  on  them  in  the  following.  A 
comparison  of  the  Helmberger-Hadley  (H-H)  and  Mueller-Murphy  (M-M)  source 
models  for  these  three  events  is  shown  in  Figure  1.  The  first  point  to  note 
is  that  even  though  the  two  formalisms  are  mathematically  different,  they  have 
very  similar  corner  frequencies  and  spectral  slopes  near  1  Hz.  The  similarity 
of  these  features  extends  over  a  wide  range  of  yield,  from  155  Kt  (SCOTCH)  to 
1300  Kt  (BOXCAR).  This  similarity  suggests  that  Mueller  and  Murphy  (1971)  and 
Barker  et  al.  (19^5)  were  resolving  the  same  yield  scaling  information  and 
simply  modeling  it  with  different  formalisms.  It  appears  that  the  best 
agreement  between  the  two  source  models  would  be  obtained  if  the 
Helmberger-Hadley  model  were  given  less  overshoot.  The  Mueller-Murphy  source 
models  have  larger  spectral  displacement  amplitudes  than  the  Helmberger-Hadley 
source  models  for  all  three  events.  As  we  show  in  the  following,  this  means 
that  utilization  of  the  Mueller-Murpby  source  requires  slightly  higher  values 
of  t*.  It  is  also  important  to  note  that,  even  though  the  Helmberger-Hadley 
source  ultimately  has  a  higher  spectral  falloff  rate,  it  does  not  become 
apparent  until  frequencies  exceeding  Hz.  Thus,  this  difference  between  the 
sources  is  not  significant  for  modeling  data  band-limited  to  frequencies  lower 
than  this. 

The  differences  between  the  two  source  formalisms  are  further 
illustrated  in  Figure  2.  The  frequency  dependence  of  the  spectral  slope  is 
displayed  at  the  top.  The  two  source  models  have  almost  identical  slopes 


Figure  1.  Predicted  far-field  source  displacement  spectra  for  BOXCAR, 

HALFBEAK,  and  SCOTCH.  The  Mueller -Mur phy  (1971)  source  models  are  based 
on  the  announced  yields.  The  Helmberger-Hadley  (1981)  source  models  are 
based  on  the  near-field  modeling  results  of  Barker  et  al.  (1985). 
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Figure  2.  Top:  spectral  slope  for  the  displacement  spectra  shown  in  Figure 
1.  Bottom:  the  spectral  ratio  of  the  Helraberger-Hadley  to 
Mueller-Murphy  source  spectra  in  Figure  1. 


between  0.5  and  3  Hz.  The  spectral  slope  of  the  Helmberger-Hadley  model  does 
not  drop  below  that  of  the  Mueller-Murphy  source  until  frequencies  from  2  to  4 
Hz.  The  slope  of  the  Helmberger-Hadley  spectrum  does  not  reach  -3  until  10 
Hz.  This  illustrates  the  danger  of  attempting  to  correct  spectra  for  source 
effects  hy  simply  multiplying  by  fn.  The  spectral  ratios  of  the  source 
functions  are  shown  at  the  bottom  of  Figure  2.  They  illustrate  again  that  the 
Mueller-Murphy  model  consistently  predicts  slightly  higher  amplitudes  at  all 
frequencies.  The  spectral  ratios  also  demonstrate  that  the  ratio  of  the 
source  models  vary  with  frequency  by  no  more  than  a  factor  of  3,  even  up  to  10 
Hz.  Furthermore,  the  ratio  of  the  source  models  for  HALFBEAK  and  BOXCAR 
varies  by  less  than  a  factor  of  2  between  0.5  and  4  Hz,  thus  making  these  two 
source  models  difficult  to  distinguish  from  one  another  using  observed  data  in 
this  band. 

In  the  spectral  shape  and  decay  analysis,  absolute  amplitude  information 
is  not  utilized,  so  the  baseline  difference  between  the  two  source  models  has 
no  effect.  It  is  apparent  in  Figure  2  that  given  spectral  data  in  the  0.5-4 
Hz  band,  the  difference  between  correcting  the  observed  spectrum  for  the 
Mueller-Murphy  source  model  versus  the  Helmberger-Hadley  model  will  depend 
somewhat  on  the  yield  of  the  event,  being  negligible  for  large  events  and 
increasing  slightly  for  small  events.  Also,  for  small  events,  the 
Helmberger-Hadley  source  correction  will  actually  lead  to  lower  Q  values.  Just 
the  opposite  effect  of  simply  correcting  for  f”^  instead  of  f-^. 


RESULTS 


Time-Domain  Amplitudes.  We  begin  this  modeling  study  by  considering  the 

amplitude  and  magnitude  data  of  Lay  (1985).  Since  we  do  not  have  accurate 

spectra  for  the  WWSSN  data,  we  model  it  in  the  time  domain.  We  compute 

synthetic  seismograms  at  a  distance  of  50°  using  the  near-field  source  models 

convolved  with  the  appropriate  attenuation  model  and  the  WWSSN  short-period 

instrument  response.  The  synthetic  seismograms  include  the  direct  P  and  pP 

arrivals.  The  pP  delay  times  and  amplitudes  were  measured  using  a  complete 

waveform  analysis  called  intercorrelation  (Lay,  1985).  Figure  3  illustrates 

the  procedure  for  determining  the  attenuation  models  from  the  amplitude  data 

for  the  event  HALFBEAK.  In  the  case  shown,  we  found  values  for  t*  in  the 

Futterraan  operator  (top)  and  t  in  the  ?4inster  operator  that  match  the 

m 

observed  average  amplitude.  We  consider  three  different  values  of  t*  in  the 
Minster  model  (1.0,  1.15  and  1.3  sec)  as  illustrated  in  the  bottom  three  rows 
of  the  figure. 

The  synthetic  waveforms  for  the  Helmberger-Hadley  source  are  different 
enough  to  allow  us  to  select  a  preferred  value  of  t*.  For  comparison,  shown 
at  the  right  of  Figure  3  are  representative  observed  waveforms  for  the  event 
HALFBEAK.  For  the  Helmberger-Hadley  source  with  the  Minster  Q  model,  the 
first  synthetic  waveform  (t*=l  sec)  resembles  the  observations  most  closely. 
The  period  of  the  first  upswing  is  too  short  and  its  relative  amplitude  too 
high  in  the  lower  two  synthetic  waveforms.  Also,  the  shoulder  in  the  second 
upswing,  which  is  associated  with  the  pP  arrival,  is  not  clear  in  the  data. 
There  are  trade-offs  here  with  the  pP  parameters,  but  we  are  assuming  that  the 
values  measured  by  Lay  (1985)  are  correct.  For  the  Mueller-Murphy  source,  it 
is  clear  that  there  is  a  strong  trade-off  between  the  two  parameters  in  the 


Figure  3.  Teleseismic  WWSSN  short-period  synthetic  seismograms  for  HALFBEAK. 
The  source  functions  are  shown  in  Figure  1.  The  attenuation  models  are 
those  that  produce  synthetic  amplitudes  which  match  the  observed 
amplitudes.  The  synthetic  seismograms  include  a  pP  arrival  (given  by 
Lay,  1985).  Also  given  for  each  synthetic  seismogram  is  the  synthetic 
magnitude.  Shown  at  the  right  are  representative  WWSSN  short-period 
teleseismic  observed  waveforms  for  HALFBEAK. 


P.MKM  ■  -M.  F  to  n  M  !-y  »  to  ^  to - to  ^4*  -  ■  r  a  ''•to  r  -to  -  -  •■  ■  - 

Minster  Q  operator.  There  are  no  substantial  differences  between  the 

* 

predicted  waveforms.  Fortunately,  a  preferred  value  of  t  can  be  established 
through  consideration  of  the  magnitude  combined  with  the  amplitude 
information. 

The  analysis  of  the  magnitude  data  of  Lay  (1995)  is  carried  out  in  the 
same  way  as  that  of  the  amplitude  values.  We  determine  values  for  the 
parameters  in  the  attenuation  models  that  match  given  by 

(1)  =  log (  A/T  )  +  P(A) 

where  A  is  amplitude  in  millimicrons,  T  is  period,  and  P(A)  is  the  distance 
correction  given  by  Veith  and  Clawson  (1972). 

Modeling  studies  like  the  one  illustrated  in  Figure  3  were  carried  out 
on  the  five  events  with  announced  yield  using  the  Mueller-Murphy  source  and  on 
the  six  events  with  near-field  models  using  the  Helmberger-Hadley  source.  The 
results  were  averaged  and  are  presented  in  Table  2.  The  average  Futterman  t* 
values  for  the  Helmberger-Hadley  source  are  0.79  and  0.75  sec  based  on  the 
amplitude  and  magnitude  observations,  respectively.  The  corresponding  values 
for  the  Mueller-Murphy  source  are  0.99  and  1.05.  The  error  bounds  given  in 
the  table  represent  one  standard  deviation.  The  result  that  the  t*  estimates 
are  about  0.2  sec  larger  for  the  Mueller-Murphy  source  compared  to  the 
Helmberger-Hadley  source  model  is  related  to  the  fact  that  the  Mueller-Murphy 
source  model  contains  significantly  more  1  Hz  energy  than  the 

Helmberger-Hadley  source  model  for  events  with  similar  yield  (see  Figures  1 

£ 

and  2).  The  uncertainties  in  the  average  t  measurements  are  generally 
smaller  for  the  Mueller-Murphy  than  for  the  Helmberger-Hadley  source  model. 
Also,  the  uncertainties  are  generally  smaller  for  the  amplitude  observations 
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Table  2 


AVERAGE  ATTENUATION  MODELS 


Futterman  Attenuation  Models: 


Murphy  Source _ Helmberger  Source 

O.99+0.OU  (M*F)  0.79+0.06  (H»F) 


t  from  Log  A  0.99+0.01*  (M»F) 


t  from 


1.05+0.08 


0.75+0.08 


Minster  Attenuation  Models: 


Murphy  Sou rc e 


tQ=1.25  seconds 


Helmberger  Source 


tQ=1.00  seconds 


I  from  Log  A  0.051+0.008  (M*M)  O.OU+O.Oll  ( H*M ) 


Tm  ^rora  '"b  0.052+0.016 


0.060+0.016 
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than  for  the  magnitude  observations. 

n 

It  is  still  necessary  to  select  a  preferred  value  of  t  appropriate  for 
Mueller-Murphy  source  model.  We  found  that  on  average  over  all  the  events 

ft 

studied  the  value  of  tQ=1.25  sec  made  it  easiest  to  simultaneously  fit  the 

amplitude  and  magnitude  observations  with  the  Mueller-Murphy  source.  We 

therefore  selected  it  as  our  preferred  value.  The  average  j  that  matched  the 

m 

log  amplitude  data  (0.051)  is  almost  identical  to  the  value  that  matched  the 

magnitude  data  (0.052).  This  indicates  that  this  source-attenuation  model 

predicts  both  the  observed  amplitude  and  period  of  the  teleseismic  waveforms. 

The  average  values  and  standard  deviations  of  t  are  given  in  Table  2.  The  t 

m  m 

measurements  (with  t*=1.0  sec)  assuming  the  Helmberger-Hadley  source  model  are 

different  for  the  log  amplitude  (0.0UU)  and  magnitude  (0.060)  data.  The 

result  that  the  magnitude  observations  give  a  larger  t  than  the  amplitude 

m 

data  indicates  that,  for  the  Minster  Q  model  that  predicts  the  observed 
amplitudes,  the  period  is  too  low  (i.e.  the  Helmberger-Hadley  waveforms  are 
too  high  in  relative  frequency  content).  We  choose  t*=1.0  sec  as  the 
preferred  value  for  the  Helmberger-Hadley  source  model  because  that  value 

ft 

represents  a  lower  bound  for  tQ.  Larger  values  will  only  result  in  synthetic 

waveforms  that  are  higher  frequency  since  a  larger  t*  will  require  a  larger  x 

0  m 

to  match  the  observed  amplitudes. 

The  Futterman  t  measurements  also  indicate  that  the  Helmberger-Hadley 
synthetic  waveforms  are  higher  in  relative  frequency  content  than  the 
Mueller-Murphy  waveforms  in  the  0.5-2  Hz  frequency  band.  (This  result  is  most 
clearly  shown  in  Figures  1  and  2  and  will  be  further  illustrated  in  the  next 

ft 

section.)  If  the  measured  t  from  log  amplitudes  is  larger  than  that  from  the 
magnitudes  (as  for  the  Helmberger-Hadley  source  model),  then  the  synthetic 
seismograms  are  too  high  in  relative  frequency.  For  the  Mueller-Murphy  source 


model,  the  Futterman  Q  model  gives  synthetic  waveforms  which  are  too 
long-period.  However,  the  Minster  attenuation  model,  which  does  not  attenuate 
high  frequencies  as  strongly,  matches  both  the  observed  amplitude  and  period. 

We  believe  that  the  models  based  on  the  log  amplitude  observations  are 
more  accurate  than  those  based  on  the  magnitude  observations.  This  is  because 
the  amplitude  data  involves  only  one  measurement  while  the  magnitudes  require 
two,  with  the  signal  period  being  occasionally  difficult  to  measure  in  the 
data.  Furthermore,  the  period  nKasurement  is  somewhat  influenced  by  pP  delay 
time.  Errors  in  the  pP  delay  time  may  alter  the  attenuation  models  that  fit 
the  magnitude  data. 

Synthetic  seismograms  for  the  two  source  models  and  the  average 
attenuation  models  (based  on  the  log  amplitude  observations)  for  the  BOXCAR, 
HALFBEAK,  and  SCOTCH  events  are  shown  in  Figure  U.  Also  given  are  the  log 
amplitude  and  magnitude  residuals.  These  are  the  difference  between  the 
observed  and  predicted  values  (a  positive  residual  means  that  the  predicted 
value  is  smaller  than  the  observed).  These  residuals  indicate  how  well  the 
various  average  attenuation  models  match  the  observed  amplitudes  and 
magnitudes.  The  log  amplitude  residuals  are  shown  in  Figure  5  for  each 
attenuation  and  source  model.  The  events  are  ordered  by  estimated  yield.  No 
trend  of  the  residuals  with  yield  is  apparent.  In  general,  the  Mueller-Murphy 
residuals  are  slightly  smaller  than  the  Helmberger-Hadley  residuals. 

To  illustrate  how  well  the  attenuation  models  based  on  the  log  amplitude 
data  predict  the  magnitudes,  the  magnitude  residuals  for  those  models  are 
given  in  Figure  6.  These  tend  to  be  larger  than  the  log  amplitude  residuals. 
This  is  not  a  surprising  result  since  the  average  attenuation  models  are  based 
on  the  amplitude  data  and  the  observed  magnitudes  are  subject  to  errors  in 
period  measurements.  Again,  no  trend  in  residual  with  yield  is  apparent. 
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Figure  U.  Teleseismic  WWSSN  short-period  synthetic  seismograms  for  BOXCAR, 
HALFBEAK,  and  SCOTCH.  The  source  functions  are  those  shown  in  Figure  1 
The  attenuation  models  are  the  average  models  that  predict  the  log 
amplitude  observations  (given  in  Table  2).  Shown  with  the  synthetic 
seismograms  are  the  log  amplitude  and  magnitude  residxials. 
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Figure  5.  Log  amplitude  residuals  for  each  of  the  attenuation-source  models. 
The  average  attenuation  models  are  from  log  amplitude  observations  and 
are  given  in  Table  2.  The  source  models  are  based  on  the  announced 
yields  (Mueller-Murpty"  source)  or  on  near-field  modeling 
(Helmberger-Hadley  source). 
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Figure  6.  Magnitude  residuals  for  each  of  the  attenuation-source  models. 

average  attenuation  models  are  from  log  amplitude  observations  and  are 
given  in  Table  2. 
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This  type  of  analysis  gives  us  an  indication  of  the  deviation  between 
the  observed  amplitude  measurements  and  the  predicted  values  based  on  the  true 
yield  and  the  average  attenuation  models.  The  residuals  for  each 
source-attenuation  model  are  generally  much  less  than  the  standard  deviations 
of  the  observed  amplitude  measurements  given  in  Table  1.  The  root  mean  square 
of  the  residuals  for  the  log  amplitude  residuals  are  O.Ol  (M*F),  0.05  (M*M), 
0.07  (H*F),  and  0.06  (H*M).  In  these  acronyms  for  the  source-attenuation 
models,  the  first  letter  refers  to  the  source  model  (H=Helmberger-Hadley , 
M=Mueller-Murphy ) ,  the  second  letter  refers  to  the  attenuation  model 
(F=Futterman ,  M=Minster),  and  the  *  refers  to  a  convolution.  The  RMS  residual 
for  the  magnitude  residuals  are  0.0?  (M*F),  0.07  (M*M),  0.08  (H*F),  and  0.11 
(H#M).  These  are  less  than  the  RMS  standard  deviations  of  the  observed  log 
amplitude  measurements  (0.10)  and  magnitude  measurements  (0.11).  Thus  the 
scatter  of  observed  amplitudes  and  magnitudes  about  the  predictions  of  our 
models  is  less  than  the  scatter  in  the  raw  amplitude  and  magnitude 
measurements. 

The  fundamental  conclusion  of  our  time-domain  amplitude  and  waveform 
modeling  study  is  that  from  WWSSN  data  either  of  the  two  source  ncdels  and 
either  of  the  two  attenuation  models  is  acceptable.  By  varying  the  free 
parameters  in  the  two  attenuation  models  within  reasonable  ranges,  we  can 
obtain  essentially  identical  predictions  of  the  amplitudes  and  waveforms.  The 
values  of  the  free  parameters  which  fit  the  data  (Table  2)  are  significant  new 
results. 

Spectral  Shape  Constraints.  The  amplitude  spectra  for  the  HALFBEAK 
synthetic  seismograms  from  Figure  ^  are  shown  in  Figure  7.  The  amplitude 
spectra  for  the  various  source-attenuation  models  are  very  similar  to  each 
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Figure  7.  Teleseisraic  WWSSN  short-period  displacement  spectra  for  the 
HALFBEAK  synthetic  seismograms  shown  in  Figure  U. 


other.  It  would  certainly  be  difficult  to  distinguish  them  with  band-limited 
data  like  that  from  the  WWSSN  short-period  instrument,  even  if  it  were 
digital.  The  maximum  amplitude  of  the  central  peak  and  its  initial  falloff 
rate  in  the  data  are  apparently  well  resolved  since  we  obtained  almost 
identical  results  for  all  four  models.  Close  inspection  of  Figure  7  shows 
that  the  Helmberger-Hadley  synthetic  seismograms  have  a  higher  relative 
frequency  content  than  the  Mueller-Murphy  waveforms  in  this  band.  Also,  the 
Minster  Q  model  synthetic  seismograms  are  higher  frequency  than  the  Futterraan 

Q  model  synthetic  seismograms.  Based  on  the  result  that  the  Minster  Q  model 

* 

with  t  =1.25  and  T  =0.051  for  the  Mueller-Murply  source  model  predicts  both 
the  amplitude  and  period  of  the  observed  waveforms,  we  infer  that  the  shape  of 
the  observed  spectra  nest  closely  resemble  the  corresponding  spectrum. 

However,  when  considering  the  amount  of  scatter  in  the  amplitude  observations, 
the  scatter  in  the  attenuation  measurements,  and  the  similarity  of  the  various 
attenuation-source  models;  we  feel  that  each  of  the  four  models  is 
acceptable. 

To  ultimately  select  the  best  attenuation  model,  we  must  consider 
broadband  spectral  observations  (e.g.  Der  et  al. ,  1985;  Bache  et  al. ,  1986). 
Der  &I«  (1985)  computed  synthetic  spectral  templates  and  compared  them  to 
observed  bocty  wave  spectra  of  Pahute  Mesa  events  to  find  the  attenuation  model 
that  best  matched  the  spectral  shape.  Th?y  assumed  the  von  Seggern-Blandford 
(1972)  granite  source  model.  Although  many  similar  studies  (e.g.  Der  et  al. , 
1982a;  Der  and  Lees,  1985)  indicate  that  a  frequency  dependent  attenuation 
model  is  in  general  necessary,  Der  et  al.  (1985)  found  that  the  spectral  shape 
in  the  0.5-**  Hz  band  can  be  Just  as  easily  fit  with  a  frequency  independent  Q 
model.  They  obtain  apparent  t  values  in  the  0.1*-0.5  sec  range.  In  general, 
the  apparent  value  of  t  appears  to  be  a  stable  parameter  and  differs  only 


£ 

slightly  from  the  absolute  value  of  t  (Der  et  al. ,  1985).  However,  for  the 

n 

frequency  independent  Q  model,  the  apparent  value  of  t  is  equivalent  to  the 

*  || 

absolute  value  of  t  .  We  therefore  assume  that  a  frequency  independent  t  of 

0.U5  sec  matches  the  spectral  shape  in  the  0.5-4  Hz  band  and  apply  it  as  a 
constraint  to  our  time-domain  amplitude  models. 

Bache  et  al.  (1986)  corrected  their  Pahute  Mesa  spectral  data  for  the 
Mueller-Murphy  source  spectrum  and  then  solved  for  the  attenuation  operator 
directly.  They  find  two  frequency  dependent  attenuation  models  that  match  the 
spectral  shape  data  in  the  0.5  to  6  Hz  band.  One  is  a  Minster  Q  model  (single 

n 

absorption  band)  with  t  =1.0  sec  and  x  =0.0*4  (designated  by  QN2).  The  other 

o  m 

n 

is  a  Minster  Q  model  with  t  =1.2  sec  and  x  =0.06  convolved  with  a  Futterman  Q 

0  m 

#  , 

model  with  t  =0.1  sec  (designated  by  QN1).  In  either  study,  the  goal  was  to 
find  the  attenuation  model  which  best  explained  the  spectral  falloff  over  a 
broad  range  of  frequencies.  As  we  show,  the  spectral  shapes  of  Pahute  Mesa 
events  measured  in  the  two  studies  are  essentially  the  same.  However,  the 
absolute  spectral  level  is  quite  different. 

The  left  portion  of  Figure  3  shows  the  spectral  shapes  of  the  source 
models  (for  HALFBEAK  as  an  example)  convolved  with  the  attenuation  models 
presented  in  this  study.  The  spectra  have  been  multiplied  by  frequency 
squared.  The  four  source-attenuation  models  described  in  this  paper 
necessarily  have  the  correct  spectral  amplitude  at  around  1  Hz  since  they  are 
based  on  matching  the  time-domain  amplitudes  on  the  WWSSN  short-period 
instrument.  Shown  as  dotted  curves  are  models  which  fit  the  spectral  shape  in 
this  band  from  Der  et  al.  (I9B5)  and  Bache  et  al.  (1986).  The  curve  from  Der 
et  al.  (1985)  is  a  von  Seggern-Blandford  granite  source  with  a  frequency 
independent  t*  of  0.45  sec.  The  curve  from  Bache  et  al.  (1986)  is  a 

Muel ler-Murphy  source  with  the  attenuation  model  QN1.  All  models  are  plotted 
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Figure  8.  Left  side:  Comparison  of  the  amplitude  spectrum  of  HALFBEAK  source 
model  convolved  with  the  attenuation  model  of  the  t*  source-attenuation 
models  presented  in  this  paper  with  the  results  of  Der  et  al»  (1985)  and 
of  Bache  et  al.  (19^6).  Right  side:  Comparison  of  the  attenuation  model 
spectrum  of  the  results  of  this  study  (only  Minster  Q  models)  with  the 
results  of  Der  et  al.  (1985)  and  of  Bache  et  al.  (1986). 
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at  their  true  absolute  levels.  We  find  that  the  results  of  Der  et  al.  (1985) 
and  Bache  et  al.  (1986)  have  similar  spectral  shapes,  but  have  different 
absolute  levels.  If  we  use  the  spectral  shape  data  as  a  constraint  for  our 
source-attenuation  models,  we  find  that  either  source  model  convolved  with  the 
appropriate  Minster  Q  model  satisfies  the  spectral  shape  as  well  as  the 
amplitude  data.  The  slope  of  the  source-attenuation  models  between  1  and  U  Hz 
as  determined  from  linear  regressions  are  similar  for  the  four  models;  0.U7 
for  QN1,  O.56  for  Der  et  al  (1985),  0.U4  for  M»M,  and  0.1*1  for  H*M.  Each  of 
these  fits  have  regression  coefficients  exceeding  0.98.  The  differences  in 
slope  is  partially  due  to  differences  in  the  source  spectral  slope  of  the  von 
Seggern-Blandford,  Mueller-Murply ,  and  Helmberger-Hadley  source  irodels. 

The  right  portion  of  Figure  8  shows  the  attenuation  model  spectrum  for 
the  results  of  this  study  (only  the  Minster  Q  models)  compared  with  the 
results  of  Der  et  al.  (1985)  and  of  Bache  et  al.  (1086).  Each  of  the  five 
models  presented  here  have  similar  shapes  over  the  0.5  to  U  Hz  band,  thus  each 
model  is  considered  to  satisfy  the  spectral  shape  data.  It  is  surprising  that 
the  attenuation  models  of  this  study  and  of  Bache  et  al.  (1986)  are  so  similar 
considering  that  the  data  analyzed  was  in  such  different  forms.  This  is 
somewhat  fortuitous  since  our  amplitude  observations  provide  no  resolution  at 
higher  frequencies.  However,  the  inclusion  of  the  spectral  shape  data  as  a 
constraint  to  our  analysis  confirms  that  our  determination  of  the  frequency 
dependence  of  Q  is  correct.  Because  of  the  different  approaches  taken,  it  is 
convincing  that  these  attenuation  models  are  appropriate  for  average  Pahute 
Mesa  observations  in  this  band. 


YIELD  ESTIMATION 
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If  we  can  assume  that  the  average  attenuation  models  given  here  are 
accurate  and  that  the  observed  log  amplitudes  are  accurate,  we  can  estimate 
the  yield  of  this  set  of  explosions.  The  yields  for  the  eight  events  were 
estimated  from  each  of  the  four  source-attenuation  models  from  the  log 
amplitude  data  and  are  given  in  Table  3*  To  estimate  the  yield  using  the 
Helmberger-Hadley  source  rrcdel  requires  that  we  define  yield  scaling  relations 
for  values  of  K  and  These  were  based  on  theoretical  and  empirical  results 

outlined  by  Barker  et  al.  (1985).  These  are 


(2)  K  =  C1  Y ~T/36^  and 


(3)  r  =  c0  Y0*91, 


where  Y  is  yield  and  and  Cg  are  constants.  The  constants  are  obtained  by 
substituting  the  announced  yield  and  near-field  values  for  BOXCAR  into 
Equations  (2)  and  (3).  Hence,  C-^  is  26.21  and  Cg  is  1.76x10^. 

The  log  yield  residuals  (estimated  minus  announced)  for  the  five  events 
with  announced  yields  are  given  in  Figure  9»  Combining  the  log  amplitude 
residuals  shown  in  Figure  5  with  the  log  yield  residuals  in  Figure  9,  we  can 
estimate  the  relation  between  the  error  in  estimated  yield  and  the  error  in 
observed  amplitude.  For  the  Mueller-Murphy  source  model  with  the  Minster 
attenuation  model, 


(U)  Alog  Y  =  -1.3  Alog  A. 
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Table  3 


YIELD  ESTIMATES  (Kt)  FROM  THE  OBSERVED  LOO  AMPLITUDES 


Announced 


Event 

Yield 

ALMENDRO 

BENHAM 

1150 

BOXCAR 

1300 

GREELEY 

870 

HALFBEAK 

INLET 

MAST 

365 

SCOTCH 

155 

M-M 

M-M 

Futterman 

Minster 

769 

778 

1158 

1218 

1068 

1100 

910 

929 

U36 

452 

352 

349 

445 

14(49 

l4l 

134 

H-H 

H-H 

Futterman 

Minster 

655 

670 

973 

1013 

891 

922 

767 

789 

U03 

407 

325 

324 

4014 

4o8 

143 

139 

Source  Model 
Q  Model 


Residual 


X 


Figure  9*  Log  yield  residuals  for  the  five  events  with  announced  yields  for 
each  of  the  source-attenuation  models.  The  log  yield  residuals  are  the 
announced  subtracted  from  the  predicted. 


With  the  exception  of  SCOTCH,  the  yield  residuals  tend  to  become  more  positive 
with  increasing  yield  for  each  of  the  source-attenuation  models.  However, 
these  residuals  are  quite  small  and  may  be  insignificant  considering  the 
amount  of  scatter  in  the  observations  and  the  attenuation  measurements. 

Using  the  estimated  yields  from  Table  3,  we  computed  the  predicted 
source  spectra  for  these  events.  Figure  10  shows  the  comparison  between  the 
estimated  source  spectra  using  the  Minster  Q  model  with  the  Mueller-Murphy 
source  and  the  predicted  source  spectra  from  the  announced  yield.  The 
difference  between  the  announced  and  the  estimated  spectra  are  shown  to  be 
minimal.  It  would  certainly  be  difficult  to  distinguish  between  the  predicted 
and  estimated  source  spectra  from  observations.  For  the  other 
source-attenuation  models,  the  differences  between  the  predicted  and  estimated 
source  spectra  are  on  the  same  level  as  those  in  Figure  10. 


DISCUSSION 


ft 

There  has  been  much  recent  discussion  as  to  the  value  of  t  at 
shorter-periods.  The  typical  value  of  t  =1.0  sec  at  long  periods  results  in 
extremely  large  attenuation  at  higher  frequencies  for  teleseismic 
observations.  Der  et  al.  (1982a)  state  that  U  Hz  amplitudes  would  be  reduced 
by  four  orders  of  magnitude  relative  to  1  Hz  amplitudes  as  the  result  of  this 

ft 

t  value.  The  results  presented  in  the  left  portion  of  Figure  8  illustrate 
this  effect.  However,  numerous  investigations  show  that  there  are  substantial 

ft 

amounts  of  U-10  Hz  energy,  far  greater  than  predicted  by  a  t  of  1  sec. 

ft 

Because  of  this  observation,  t  is  thought  to  be  much  smaller  at  higher 


frf  quencies. 


t  Amplitude, 


Using  the  preferred  source-attenuation  models  presented  in  Figure  %  we 
can  compute  the  value  of  t*  as  a  function  of  frequency.  These  results  are 

ft 

presented  in  Figure  11.  Shown  are  the  values  of  t  for  the  h 
source-attenuation  models  presented  in  this  study,  along  with  the  results  of 
Bache  et  al.  (1986)  and  Der  and  Lees  (1985).  The  Futterman  Q  models  from  this 
study  (M*F  and  H*F)  are  only  conside  ;d  appropriate  near  1  Hz.  The 
attenuation  models  M*M,  H*M,  QN1,  and  QN2  are  considered  to  match  both  the 

amplitude  and  spectral  shape  data  over  this  band.  We  see  from  this  that  the 

» 

value  of  t  around  1  Hz  is  between  0.8  and  1.0  sec,  depending  on  the  choice  of 
source-attenuation  model.  Der  and  Lees  (1985)  proposed  a  frequency  dependent 
attenuation  model  for  the  western  United  States  that  satisfied  spectral  shape 
date  and  time-domain  waveform  constraints.  The  events  they  studied  were  deep 
earthquakes,  so  the  data  was  not  subject  to  upper  mantle  attenuation  on  the 
source  side  of  the  raypath.  Their  preferred  model  (QP  S-T)  is  shown  *'or 
comparison  in  Figure  11.  (S-T  refers  to  shield-to-tectonic  paths.)  Although 
Der  and  Lees  (1985)  argue  that  previous  work  (e.g.  Der  et  al. ,  1982b) 

ft 

demonstrate  that  the  apparent  t  from  deep  earthquakes  is  similar  to  that 
obtained  from  shallow  events  in  shield  regions,  it  is  reasonable  to  expect 
that  the  QP  S-T  curve  represents  a  lower  bound.  If  we  were  to  assume  some  t* 
contribution  from  the  source  side  of  the  raypath,  then  the  QP  S-T  curve  would 
become  consistent  with  the  other  attenuation  models. 

ft 

t  values  in  the  range  of  0. 8-1.0  sec  near  1  Hz  are  quite  consistent 
with  other  estimates  of  t*.  Burdick  et  al.  (198U)  obtain  a  value  of  0.9  sec 
from  an  analogous  near-field  to  far-field  amplitude  comparison  from  Amchitka 
nuclear  explosions.  Burdick  (1978)  conducted  a  source  canceling  experiment 
comparing  sP  with  sS  phases  and  obtained  a  value  of  t*=1.3  sec  at  about  0.5  Hz 
for  an  event  in  Southern  California  (presumably  similar  to  NTS).  This  value 


is  consistent  with  the  M*M  and  QN1  models.  Burdick  (1985)  estimated  a  global 


average  t  of  about  1  sec  near  1  Hz  from  a  comparison  of  ScP  and  ScS  phases. 

* 

Finally,  Burdick  and  Grand  (1985)  obtain  t  =0.99  sec  at  around  1  Hz  from  the 
comparison  of  sP  and  sS  waves  from  Asia. 

To  illustrate  the  effect  of  the  absolute  value  of  t  on  the  observed 
amplitudes,  we  computed  synthetic  seismograms  for  the  Mueller-Murphy  source 

model  with  various  Q  models  with  similar  spectral  shapes  in  the  0.5-4  Hz  band. 

* 

We  assumed  the  M»M  source-attenuation  model  (this  paper),  along  with  tQ=0.75 

* 

sec  and  1^=0.025,  and  a  frequency  independent  t  =0.42  sec  (similar  to  Der  et 
al .  ,  1985).  Each  of  these  Q  models  have  the  same  spectral  slope  (0.57)  in  the 
0.5-4  Hz  band  (linear  regression  coefficients  exceed  0.97).  The  results  for 
the  HALFBEAK  event  are  presented  in  Figure  12.  The  predicted  amplitude  using 
t*=0 .42  sec  is  nearly  a  factor  of  5  larger  than  observed,  while  the  predicted 
amplitude  for  t*=1.25  and  1^=0.051  is  similar  to  the  observed  value. 
Furthermore,  the  synthetic  seismogram  with  t  =0.42  sec  is  too  high  in  relative 
frequency  content  compared  to  representative  observations  (shown  at  the  right 
in  Figure  3).  This  example  illustrates  that  the  use  of  spectral  shape  data 
without  absolute  amplitude  information  can  lead  to  incorrect  Q  models.  We 
have  previously  shown  that  the  use  of  absolute  amplitude  information  without 
spectral  shape  data  can  also  lead  to  Incorrect  results.  For  example,  the  M*M 
and  M*F  models  have  the  same  absolute  amplitude  level  near  1  Hz,  but  the 
spectral  slopes  in  the  0.5-4  Hz  band  are  very  different.  It  is  thus  necessary 
to  use  both  absolute  amplitude  information  and  spectral  shape  or  decay 
information  to  obtain  reliable  attenuation  models. 

A  direct  application  of  our  results  is  that  we  are  able  to  predict  the 
relation  between  yield  and  magnitude  based  upon  the  appropriate  attenuation 
models  for  a  given  site.  We  do  so  here  assuming  the  Mueller-Murphy  source  and 
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Figure  12.  Teleselsmic  WWSSN  short-period  synthetic  seismograms  for  RALFBEAK. 
Shown  are  the  synthetic  waveforms  Tor  the  Mueller-Murphy  source  model 
with  the  Minster  attenuation  model  presented  in  this  paper  compared  to 
other  attenuation  model t  with  similar  spectral  slopes  but  with  different 
absolute  levels. 
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the  Minster  Q  model  (M*M) .  Synthetic  seismograms  were  computed  as  in  the 
above  analysis  for  yields  ranging  from  50  to  1300  Kt.  The  pP  delay  times  were 
determined  from  the  average  pP  delay  -  depth  relation  determined  from  the 
results  of  Lay  (1985)  and  assuming  the  predicted  depth  -  yield  dependence  from 
Mueller  and  Murphy  (1971).  The  predicted  magnitudes  are  given  in  Figure  13. 
The  linear  relation  (with  regression  coefficient  of  1 )  is  given  by 

(5)  =  3.695  +  0.8019  Log  Y. 

For  the  predicted  amplitudes,  the  relation  is  given  by 

(6)  Log  A  =  0.3862  +  0.7252  Log  Y. 

The  difference  in  slopes  in  Equations  5  and  6  is  due  to  the  period  correction 
included  in  the  magnitude  calculation.  Equation  5  is  almost  identical  to  that 
obtained  by  Lay  (1985)  from  a  regression  of  observed  magnitudes  for  the  five 
events  with  announced  yields.  The  observed  magnitudes  are  also  given  in 
Figure  13  and  are  in  good  agreement  with  the  predicted  m^-yield  relation. 

We  have  not  yet  addressed  the  possibility  of  constraining  the  source 
models  with  longer-period  observations.  As  apparent  in  Figure  1,  very  precise 
constraints  on  10  sec  spectral  levels  could  serve  to  discriminate  between  the 
source  models.  However,  the  long-period  levels  for  the  Helmberger-Hadley 
source  could  be  increased  simply  by  decreasing  the  overshoot  parameter  B, 
which  was  arbitrarily  set  to  1  in  the  near-field  modeling.  The  near-field 
waveforms  are  essentially  unaffected  if  the  value  of  B  is  set  to  0.5,  which 
gives  ^  a  factor  of  1.7  larger.  Smaller  values  of  B  will  increase  4* 
further,  but  will  begin  to  affect  the  spectrum  in  the  near-field  band  (1-5 
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Figure  13.  The  predicted  n^-yield  relation  for  WWSSN  short-period 

determinations  for  the  Mueller-Murplv  source  and  Minster  Q  model.  Sh 
are  the  observed  magnitudes  for  the  five  events  with  announced  yields 
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Hz),  and  thus  change  the  waveform  matches.  For  now  we  will  limit  our 


discussion  to  B>0.5,  but  acknowledge  the  possibility  of  even  smaller  values  of 
B.  Setting  the  value  of  B  at  0.5  will  increase  the  WWSSN  short-period  (0.5-2 
Hz)  amplitudes  by  less  than  10%,  thus  our  high  frequency  attenuation  operator 
modeling  will  be  relatively  unchanged  by  any  long-period  modifications. 

Several  studies  have  obtained  estimates  of  long-period  source  spectral 


levels  for  Pahute  Mesa  tests  from  analysis  of  surface  waves  (Stevens,  1985; 
Given  and  Mellman,  1985).  These  estimates  are  compared  with  the  predictions 
of  the  Mueller-Murphy  source  model,  as  well  as  the  near-field  modeling  results 
of  Barker  et  al.  (1985)  with  B=1 .0  in  Table  4.  On  average,  the  surface  wave 


results  are  consistent  with  the  Mueller-Murphy  source  predictions,  but  are  a 
factor  of  3-4  larger  than  the  Helmberger-Hadley  source  predictions.  Allowing 
for  a  factor  of  1.7  increase  in  <P  by  decreasing  B  in  the  near-field  models 
still  leaves  a  factor  of  about  2  discrepancy  between  the  surface  wave 
estimates  and  the  results  of  Barker  et  al .  (1985).  While  this  favors  the 
Mueller-Murphy  source  model,  further  work  is  needed  to  establish  the 
confidence  levels  on  the  various  estimates  before  a  definitive  conclusion  can 
be  drawn. 
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Table  4 


ESTIMATED  EXPLOSION  SOURCE  STRENGTH  ¥  (xlO10  cm3) 

00 


Given  and 

Stevens , 

M-M  * 

H-H  ( B=1 )  * 

Event 

Mellman,  1985 

1986 

Ann.  Yield 

Barker  et  al,  1985 

ALMENDRO 

17.8 

19.0 

(20.0) 

5.7 

BENHAM 

— 

— 

26.9 

(9.6) 

BOXCAR 

— 

— 

29.5 

12.0 

GREELEY 

— 

— 

21.7 

(7.6) 

HALFBEAK 

— 

— 

11.2 

3-8 

INLET 

— 

— 

(10.9) 

3-2 

MAST 

15.5 

12.0 

(13.2) 

4.7 

SCOTCH 

8.9 

9.1 

5.9 

1.3 

*  The  numbers  in  parentheses  ( )  are  estimated  from  the  predicted 
yields  assuming  the  Minster  Q  models. 


CONCLUSIONS 


Attenuation  models  for  short-period  P  waves  from  Pahute  Mesa  nuclear 

explosions  have  been  developed  by  time-domain  modeling  of  teleseisraic 

waveforms.  Trade-offs  between  source  and  attenuation  models  were  explored 

using  two  explosion  source  representations  and  two  attenuation  models,  one  of 

which  includes  frequency  dependence.  For  the  Mueller-Murphy  source  model, 

_  o 

which  has  an  asymptotic  high  frequency  spectral  falloff  of  f  ,  yield  scaling 

is  used  to  determine  the  source  models.  This  source  leads  to  frequency 

independent  t*  values  of  about  1.0  sec,  or  frequency  dependent  absorption  band 

models  with  t*=1.25  sec  and  t  (high  frequency  roll  off)  of  about  0.05.  The 
°  m 

Helmberger-Hadley  source,  which  has  a  high  frequency  spectral  decay  of  f“3,  is 

constrained  by  near-field  waveform  modeling.  This  source  results  in  frequency 

independent  models  with  t*  values  of  about  0.8  sec,  or  frequency  dependent 

models  with  t*=1.0  sec  and  t  of  about  0.04.  Spectral  shape  information  from 
°  m 

higher  frequency  data  is  well  matched  by  the  absorption  band  attenuation 
models  for  either  source  representation.  We  find  that  the  use  of  both 
spectral  shape  data  and  absolute  amplitude  information  is  necessary  to 
determine  reliable  attenuation  models.  Simple  correction  of  observed  spectra 
by  assuming  the  high-frequency  asymptotic  spectral  decay  rates  of  these  source 
models  in  the  short-period  band  can  lead  to  erroneous  conclusions  regarding 
the  attenuation  operator.  The  two  source  models  are  sufficiently  similar  in 
the  0.5-4  Hz  band  that  it  is  not  possible  to  select  the  most  appropriate 
model.  Long-period  measurements  tend  to  favor  the  higher  long-period  spectral 
amplitudes  predicted  by  the  Mueller-Murphy  source  model.  Average  attenuation 
models  for  each  source  model  can  be  used  to  accurately  predict  the  yields  for 


other  events 
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